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