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 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 } 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 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 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 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 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 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
