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
17U_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}
73static 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
81double 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
92double 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
104double 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
111U_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
189UBool calcSieveTime = FALSE;
190double meanSieveTime = 0.0;
191double meanSieveME = 0.0;
192
193U_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