icuSources/test/perf/howExpensiveIs/sieve.cpp

   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 }