| OLD | NEW |
| (Empty) |
| 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 *margi
nOfError) { | |
| 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 } | |
| OLD | NEW |