1 /*
2 **********************************************************************
3 * Copyright (c) 2011-2012,International Business Machines
4 * Corporation and others. All Rights Reserved.
5 **********************************************************************
6 */
7
8 #include "unicode/utimer.h"
9 #include <stdio.h>
10 #include <math.h>
11 #include <stdlib.h>
12
13 #include "sieve.h"
14
15 /* prime number sieve */
16
uprv_calcSieveTime()17 U_CAPI double uprv_calcSieveTime() {
18 #if 1
19 #define SIEVE_SIZE U_LOTS_OF_TIMES /* standardized size */
20 #else
21 #define SIEVE_SIZE <something_smaller>
22 #endif
23
24 #define SIEVE_PRINT 0
25
26 char sieve[SIEVE_SIZE];
27 UTimer a,b;
28 int i,k;
29
30 utimer_getTime(&a);
31 for(int j=0;j<SIEVE_SIZE;j++) {
32 sieve[j]=1;
33 }
34 sieve[0]=0;
35 utimer_getTime(&b);
36
37
38 #if SIEVE_PRINT
39 printf("init %d: %.9f\n", SIEVE_SIZE,utimer_getDeltaSeconds(&a,&b));
40 #endif
41
42 utimer_getTime(&a);
43 for(i=2;i<SIEVE_SIZE/2;i++) {
44 for(k=i*2;k<SIEVE_SIZE;k+=i) {
45 sieve[k]=0;
46 }
47 }
48 utimer_getTime(&b);
49 #if SIEVE_PRINT
50 printf("sieve %d: %.9f\n", SIEVE_SIZE,utimer_getDeltaSeconds(&a,&b));
51
52 if(SIEVE_PRINT>0) {
53 k=0;
54 for(i=2;i<SIEVE_SIZE && k<SIEVE_PRINT;i++) {
55 if(sieve[i]) {
56 printf("%d ", i);
57 k++;
58 }
59 }
60 puts("");
61 }
62 {
63 k=0;
64 for(i=0;i<SIEVE_SIZE;i++) {
65 if(sieve[i]) k++;
66 }
67 printf("Primes: %d\n", k);
68 }
69 #endif
70
71 return utimer_getDeltaSeconds(&a,&b);
72 }
comdoub(const void * aa,const void * bb)73 static int comdoub(const void *aa, const void *bb)
74 {
75 const double *a = (const double*)aa;
76 const double *b = (const double*)bb;
77
78 return (*a==*b)?0:((*a<*b)?-1:1);
79 }
80
midpoint(double * times,double i,int n)81 double midpoint(double *times, double i, int n) {
82 double fl = floor(i);
83 double ce = ceil(i);
84 if(ce>=n) ce=n;
85 if(fl==ce) {
86 return times[(int)fl];
87 } else {
88 return (times[(int)fl]+times[(int)ce])/2;
89 }
90 }
91
medianof(double * times,int n,int type)92 double medianof(double *times, int n, int type) {
93 switch(type) {
94 case 1:
95 return midpoint(times,n/4,n);
96 case 2:
97 return midpoint(times,n/2,n);
98 case 3:
99 return midpoint(times,(n/2)+(n/4),n);
100 }
101 return -1;
102 }
103
qs(double * times,int n,double * q1,double * q2,double * q3)104 double qs(double *times, int n, double *q1, double *q2, double *q3) {
105 *q1 = medianof(times,n,1);
106 *q2 = medianof(times,n,2);
107 *q3 = medianof(times,n,3);
108 return *q3-*q1;
109 }
110
uprv_getMeanTime(double * times,uint32_t * timeCount,double * marginOfError)111 U_CAPI double uprv_getMeanTime(double *times, uint32_t *timeCount, double *marginOfError) {
112 double q1,q2,q3;
113 int n = *timeCount;
114
115 /* calculate medians */
116 qsort(times,n,sizeof(times[0]),comdoub);
117 double iqr = qs(times,n,&q1,&q2,&q3);
118 double rangeMin= (q1-(1.5*iqr));
119 double rangeMax = (q3+(1.5*iqr));
120
121 /* Throw out outliers */
122 int newN = n;
123 #if U_DEBUG
124 printf("iqr: %.9f, q1=%.9f, q2=%.9f, q3=%.9f, max=%.9f, n=%d\n", iqr,q1,q2,q3,(double)-1, n);
125 #endif
126 for(int i=0;i<newN;i++) {
127 if(times[i]<rangeMin || times[i]>rangeMax) {
128 #if U_DEBUG
129 printf("Removing outlier: %.9f outside [%.9f:%.9f]\n", times[i], rangeMin, rangeMax);
130 #endif
131 times[i--] = times[--newN]; // bring down a new value
132 }
133 }
134
135 #if U_DEBUG
136 UBool didRemove = false;
137 #endif
138 /* if we removed any outliers, recalculate iqr */
139 if(newN<n) {
140 #if U_DEBUG
141 didRemove = true;
142 printf("removed %d outlier(s), recalculating IQR..\n", n-newN);
143 #endif
144 n = newN;
145 *timeCount = n;
146
147 qsort(times,n,sizeof(times[0]),comdoub);
148 double iqr = qs(times,n,&q1,&q2,&q3);
149 rangeMin= (q1-(1.5*iqr));
150 rangeMax = (q3+(1.5*iqr));
151 }
152
153 /* calculate min/max and mean */
154 double minTime = times[0];
155 double maxTime = times[0];
156 double meanTime = times[0];
157 for(int i=1;i<n;i++) {
158 if(minTime>times[i]) minTime=times[i];
159 if(maxTime<times[i]) maxTime=times[i];
160 meanTime+=times[i];
161 }
162 meanTime /= n;
163
164 /* caculate standard deviation */
165 double sd = 0;
166 for(int i=0;i<n;i++) {
167 #if U_DEBUG
168 if(didRemove) {
169 printf("recalc %d/%d: %.9f\n", i, n, times[i]);
170 }
171 #endif
172 sd += (times[i]-meanTime)*(times[i]-meanTime);
173 }
174 sd = sqrt(sd/((double)n-1.0));
175
176 #if U_DEBUG
177 printf("sd: %.9f, mean: %.9f\n", sd, meanTime);
178 printf("min: %.9f, q1=%.9f, q2=%.9f, q3=%.9f, max=%.9f, n=%d\n", minTime,q1,q2,q3,maxTime, n);
179 printf("iqr/sd = %.9f\n", iqr/sd);
180 #endif
181
182 /* 1.960 = z sub 0.025 */
183 *marginOfError = 1.960 * (sd/sqrt((double)n));
184 /*printf("Margin of Error = %.4f (95%% confidence)\n", me);*/
185
186 return meanTime;
187 }
188
189 UBool calcSieveTime = FALSE;
190 double meanSieveTime = 0.0;
191 double meanSieveME = 0.0;
192
uprv_getSieveTime(double * marginOfError)193 U_CAPI double uprv_getSieveTime(double *marginOfError) {
194 if(calcSieveTime==FALSE) {
195 #define SAMPLES 50
196 uint32_t samples = SAMPLES;
197 double times[SAMPLES];
198
199 for(int i=0;i<SAMPLES;i++) {
200 times[i] = uprv_calcSieveTime();
201 #if U_DEBUG
202 printf("sieve: %d/%d: %.9f\n", i,SAMPLES, times[i]);
203 #endif
204 }
205
206 meanSieveTime = uprv_getMeanTime(times, &samples,&meanSieveME);
207 calcSieveTime=TRUE;
208 }
209 if(marginOfError!=NULL) {
210 *marginOfError = meanSieveME;
211 }
212 return meanSieveTime;
213 }
214