icuSources/test/perf/howExpensiveIs/sieve.cpp

   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 }