17 const double pi = 3.14159265358979323846264338328;
19 void release_g_lfact_table()
22 munmap(g_lfact_table, ISOSPEC_G_FACT_TABLE_SIZE*
sizeof(
double));
28 double* alloc_lfact_table()
32 ret =
reinterpret_cast<double*
>(mmap(
nullptr,
sizeof(
double)*ISOSPEC_G_FACT_TABLE_SIZE, PROT_READ | PROT_WRITE, MAP_ANONYMOUS | MAP_PRIVATE, -1, 0));
34 ret =
reinterpret_cast<double*
>(calloc(ISOSPEC_G_FACT_TABLE_SIZE,
sizeof(
double)));
36 std::atexit(release_g_lfact_table);
40 double* g_lfact_table = alloc_lfact_table();
43 double RationalApproximation(
double t)
47 double c[] = {2.515517, 0.802853, 0.010328};
48 double d[] = {1.432788, 0.189269, 0.001308};
49 return t - ((c[2]*t + c[1])*t + c[0]) /
50 (((d[2]*t + d[1])*t + d[0])*t + 1.0);
53 double NormalCDFInverse(
double p)
57 return -RationalApproximation( sqrt(-2.0*log(p)) );
59 return RationalApproximation( sqrt(-2.0*log(1-p)) );
62 double NormalCDFInverse(
double p,
double mean,
double stdev)
64 return mean + stdev * NormalCDFInverse(p);
67 double NormalCDF(
double x,
double mean,
double stdev)
69 x = (x-mean)/stdev * 0.7071067811865476;
72 double a1 = 0.254829592;
73 double a2 = -0.284496736;
74 double a3 = 1.421413741;
75 double a4 = -1.453152027;
76 double a5 = 1.061405429;
86 double t = 1.0/(1.0 + p*x);
87 double y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x);
89 return 0.5*(1.0 + sign*y);
92 double NormalPDF(
double x,
double mean,
double stdev)
94 double two_variance = stdev * stdev * 2.0;
95 double delta = x-mean;
96 return exp( -delta*delta / two_variance ) / sqrt( two_variance * pi );