1103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius/* 264339d36f8bd4db5025fe2988eda22b491a9219cFredrik Roubert *********************************************************************** 364339d36f8bd4db5025fe2988eda22b491a9219cFredrik Roubert * Copyright (C) 2016 and later: Unicode, Inc. and others. 464339d36f8bd4db5025fe2988eda22b491a9219cFredrik Roubert * License & terms of use: http://www.unicode.org/copyright.html#License 564339d36f8bd4db5025fe2988eda22b491a9219cFredrik Roubert *********************************************************************** 664339d36f8bd4db5025fe2988eda22b491a9219cFredrik Roubert *********************************************************************** 754dcd9b6a06071f647dac967e9e267abb9410720Craig Cornelius * Copyright (c) 2011-2012,International Business Machines 8103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius * Corporation and others. All Rights Reserved. 964339d36f8bd4db5025fe2988eda22b491a9219cFredrik Roubert *********************************************************************** 10103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius */ 11103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 12103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#include "unicode/utimer.h" 13103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#include <stdio.h> 14103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#include <math.h> 15103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#include <stdlib.h> 16103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 17103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#include "sieve.h" 18103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 19103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius/* prime number sieve */ 20103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 21103e9ffba2cba345d0078eb8b8db33249f81840aCraig CorneliusU_CAPI double uprv_calcSieveTime() { 22103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#if 1 23103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#define SIEVE_SIZE U_LOTS_OF_TIMES /* standardized size */ 24103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#else 25103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#define SIEVE_SIZE <something_smaller> 26103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#endif 27103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 28103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#define SIEVE_PRINT 0 29103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 30103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius char sieve[SIEVE_SIZE]; 31103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius UTimer a,b; 32103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius int i,k; 33103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 34103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius utimer_getTime(&a); 35103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius for(int j=0;j<SIEVE_SIZE;j++) { 36103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius sieve[j]=1; 37103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius } 38103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius sieve[0]=0; 39103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius utimer_getTime(&b); 40103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 41103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 42103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#if SIEVE_PRINT 43103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius printf("init %d: %.9f\n", SIEVE_SIZE,utimer_getDeltaSeconds(&a,&b)); 44103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#endif 45103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 46103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius utimer_getTime(&a); 47103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius for(i=2;i<SIEVE_SIZE/2;i++) { 48103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius for(k=i*2;k<SIEVE_SIZE;k+=i) { 49103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius sieve[k]=0; 50103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius } 51103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius } 52103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius utimer_getTime(&b); 53103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#if SIEVE_PRINT 54103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius printf("sieve %d: %.9f\n", SIEVE_SIZE,utimer_getDeltaSeconds(&a,&b)); 55103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 56103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius if(SIEVE_PRINT>0) { 57103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius k=0; 58103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius for(i=2;i<SIEVE_SIZE && k<SIEVE_PRINT;i++) { 59103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius if(sieve[i]) { 60103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius printf("%d ", i); 61103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius k++; 62103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius } 63103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius } 64103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius puts(""); 65103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius } 66103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius { 67103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius k=0; 68103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius for(i=0;i<SIEVE_SIZE;i++) { 69103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius if(sieve[i]) k++; 70103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius } 71103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius printf("Primes: %d\n", k); 72103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius } 73103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#endif 74103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 75103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius return utimer_getDeltaSeconds(&a,&b); 76103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius} 77103e9ffba2cba345d0078eb8b8db33249f81840aCraig Corneliusstatic int comdoub(const void *aa, const void *bb) 78103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius{ 79103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius const double *a = (const double*)aa; 80103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius const double *b = (const double*)bb; 81103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 82103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius return (*a==*b)?0:((*a<*b)?-1:1); 83103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius} 84103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 85103e9ffba2cba345d0078eb8b8db33249f81840aCraig Corneliusdouble midpoint(double *times, double i, int n) { 86103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius double fl = floor(i); 87103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius double ce = ceil(i); 88103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius if(ce>=n) ce=n; 89103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius if(fl==ce) { 90103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius return times[(int)fl]; 91103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius } else { 92103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius return (times[(int)fl]+times[(int)ce])/2; 93103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius } 94103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius} 95103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 96103e9ffba2cba345d0078eb8b8db33249f81840aCraig Corneliusdouble medianof(double *times, int n, int type) { 97103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius switch(type) { 98103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius case 1: 99103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius return midpoint(times,n/4,n); 100103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius case 2: 101103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius return midpoint(times,n/2,n); 102103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius case 3: 103103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius return midpoint(times,(n/2)+(n/4),n); 104103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius } 105103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius return -1; 106103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius} 107103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 108103e9ffba2cba345d0078eb8b8db33249f81840aCraig Corneliusdouble qs(double *times, int n, double *q1, double *q2, double *q3) { 109103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius *q1 = medianof(times,n,1); 110103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius *q2 = medianof(times,n,2); 111103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius *q3 = medianof(times,n,3); 112103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius return *q3-*q1; 113103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius} 114103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 11554dcd9b6a06071f647dac967e9e267abb9410720Craig CorneliusU_CAPI double uprv_getMeanTime(double *times, uint32_t *timeCount, double *marginOfError) { 116103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius double q1,q2,q3; 11754dcd9b6a06071f647dac967e9e267abb9410720Craig Cornelius int n = *timeCount; 118103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 119103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius /* calculate medians */ 120103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius qsort(times,n,sizeof(times[0]),comdoub); 121103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius double iqr = qs(times,n,&q1,&q2,&q3); 122103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius double rangeMin= (q1-(1.5*iqr)); 123103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius double rangeMax = (q3+(1.5*iqr)); 124103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 125103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius /* Throw out outliers */ 126103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius int newN = n; 127103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#if U_DEBUG 128103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius printf("iqr: %.9f, q1=%.9f, q2=%.9f, q3=%.9f, max=%.9f, n=%d\n", iqr,q1,q2,q3,(double)-1, n); 129103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#endif 130103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius for(int i=0;i<newN;i++) { 131103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius if(times[i]<rangeMin || times[i]>rangeMax) { 132103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#if U_DEBUG 13354dcd9b6a06071f647dac967e9e267abb9410720Craig Cornelius printf("Removing outlier: %.9f outside [%.9f:%.9f]\n", times[i], rangeMin, rangeMax); 134103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#endif 135103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius times[i--] = times[--newN]; // bring down a new value 136103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius } 137103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius } 138103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 13954dcd9b6a06071f647dac967e9e267abb9410720Craig Cornelius#if U_DEBUG 14054dcd9b6a06071f647dac967e9e267abb9410720Craig Cornelius UBool didRemove = false; 14154dcd9b6a06071f647dac967e9e267abb9410720Craig Cornelius#endif 142103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius /* if we removed any outliers, recalculate iqr */ 143103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius if(newN<n) { 144103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#if U_DEBUG 14554dcd9b6a06071f647dac967e9e267abb9410720Craig Cornelius didRemove = true; 14654dcd9b6a06071f647dac967e9e267abb9410720Craig Cornelius printf("removed %d outlier(s), recalculating IQR..\n", n-newN); 147103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#endif 148103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius n = newN; 14954dcd9b6a06071f647dac967e9e267abb9410720Craig Cornelius *timeCount = n; 150103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 151103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius qsort(times,n,sizeof(times[0]),comdoub); 152103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius double iqr = qs(times,n,&q1,&q2,&q3); 153103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius rangeMin= (q1-(1.5*iqr)); 154103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius rangeMax = (q3+(1.5*iqr)); 155103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius } 156103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 157103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius /* calculate min/max and mean */ 158103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius double minTime = times[0]; 159103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius double maxTime = times[0]; 160103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius double meanTime = times[0]; 161103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius for(int i=1;i<n;i++) { 162103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius if(minTime>times[i]) minTime=times[i]; 163103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius if(maxTime<times[i]) maxTime=times[i]; 164103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius meanTime+=times[i]; 165103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius } 166103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius meanTime /= n; 167103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 168103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius /* caculate standard deviation */ 169103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius double sd = 0; 170103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius for(int i=0;i<n;i++) { 171103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#if U_DEBUG 17254dcd9b6a06071f647dac967e9e267abb9410720Craig Cornelius if(didRemove) { 17354dcd9b6a06071f647dac967e9e267abb9410720Craig Cornelius printf("recalc %d/%d: %.9f\n", i, n, times[i]); 17454dcd9b6a06071f647dac967e9e267abb9410720Craig Cornelius } 175103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#endif 176103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius sd += (times[i]-meanTime)*(times[i]-meanTime); 177103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius } 17854dcd9b6a06071f647dac967e9e267abb9410720Craig Cornelius sd = sqrt(sd/((double)n-1.0)); 179103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 180103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#if U_DEBUG 181103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius printf("sd: %.9f, mean: %.9f\n", sd, meanTime); 182103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius printf("min: %.9f, q1=%.9f, q2=%.9f, q3=%.9f, max=%.9f, n=%d\n", minTime,q1,q2,q3,maxTime, n); 183103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius printf("iqr/sd = %.9f\n", iqr/sd); 184103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#endif 185103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 186103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius /* 1.960 = z sub 0.025 */ 18754dcd9b6a06071f647dac967e9e267abb9410720Craig Cornelius *marginOfError = 1.960 * (sd/sqrt((double)n)); 188103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius /*printf("Margin of Error = %.4f (95%% confidence)\n", me);*/ 189103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 190103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius return meanTime; 191103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius} 192103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 193103e9ffba2cba345d0078eb8b8db33249f81840aCraig CorneliusUBool calcSieveTime = FALSE; 194103e9ffba2cba345d0078eb8b8db33249f81840aCraig Corneliusdouble meanSieveTime = 0.0; 195103e9ffba2cba345d0078eb8b8db33249f81840aCraig Corneliusdouble meanSieveME = 0.0; 196103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 197103e9ffba2cba345d0078eb8b8db33249f81840aCraig CorneliusU_CAPI double uprv_getSieveTime(double *marginOfError) { 198103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius if(calcSieveTime==FALSE) { 199103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#define SAMPLES 50 20054dcd9b6a06071f647dac967e9e267abb9410720Craig Cornelius uint32_t samples = SAMPLES; 201103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius double times[SAMPLES]; 202103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 203103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius for(int i=0;i<SAMPLES;i++) { 204103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius times[i] = uprv_calcSieveTime(); 205103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#if U_DEBUG 20654dcd9b6a06071f647dac967e9e267abb9410720Craig Cornelius printf("sieve: %d/%d: %.9f\n", i,SAMPLES, times[i]); 207103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius#endif 208103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius } 209103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius 21054dcd9b6a06071f647dac967e9e267abb9410720Craig Cornelius meanSieveTime = uprv_getMeanTime(times, &samples,&meanSieveME); 211103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius calcSieveTime=TRUE; 212103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius } 213103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius if(marginOfError!=NULL) { 214103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius *marginOfError = meanSieveME; 215103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius } 216103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius return meanSieveTime; 217103e9ffba2cba345d0078eb8b8db33249f81840aCraig Cornelius} 218