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