18 void release_g_lfact_table()
21 munmap(g_lfact_table, ISOSPEC_G_FACT_TABLE_SIZE*
sizeof(
double));
27 double* alloc_lfact_table()
31 ret =
reinterpret_cast<double*
>(mmap(
nullptr,
sizeof(
double)*ISOSPEC_G_FACT_TABLE_SIZE, PROT_READ | PROT_WRITE, MAP_ANONYMOUS | MAP_PRIVATE, -1, 0));
33 ret =
reinterpret_cast<double*
>(calloc(ISOSPEC_G_FACT_TABLE_SIZE,
sizeof(
double)));
35 std::atexit(release_g_lfact_table);
39 double* g_lfact_table = alloc_lfact_table();
42 double RationalApproximation(
double t)
46 double c[] = {2.515517, 0.802853, 0.010328};
47 double d[] = {1.432788, 0.189269, 0.001308};
48 return t - ((c[2]*t + c[1])*t + c[0]) /
49 (((d[2]*t + d[1])*t + d[0])*t + 1.0);
52 double NormalCDFInverse(
double p)
56 return -RationalApproximation( sqrt(-2.0*log(p)) );
58 return RationalApproximation( sqrt(-2.0*log(1-p)) );
61 double NormalCDFInverse(
double p,
double mean,
double stdev)
63 return mean + stdev * NormalCDFInverse(p);
66 double NormalCDF(
double x,
double mean,
double stdev)
68 x = (x-mean)/stdev * 0.7071067811865476;
71 double a1 = 0.254829592;
72 double a2 = -0.284496736;
73 double a3 = 1.421413741;
74 double a4 = -1.453152027;
75 double a5 = 1.061405429;
85 double t = 1.0/(1.0 + p*x);
86 double y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x);
88 return 0.5*(1.0 + sign*y);
91 double NormalPDF(
double x,
double mean,
double stdev)
93 double two_variance = stdev * stdev * 2.0;
94 double delta = x-mean;
95 return exp( -delta*delta / two_variance ) / sqrt( two_variance * pi );
98 const double sqrt_pi = 1.772453850905516027298167483341145182798;
100 double LowerIncompleteGamma2(
int a,
double x)
103 double exp_minus_x = exp(-x);
107 base = 1 - exp_minus_x;
113 base = sqrt_pi * erf(sqrt(x));
120 base = base * current_s - pow(x, current_s) * exp_minus_x;
127 double InverseLowerIncompleteGamma2(
int a,
double x)
130 double p = tgamma(a);
135 v = LowerIncompleteGamma2(a, s);
140 }
while((p-l)*1000.0>p);