23 void release_g_lfact_table()
26 munmap(g_lfact_table, ISOSPEC_G_FACT_TABLE_SIZE*
sizeof(
double));
32 double* alloc_lfact_table()
36 ret =
reinterpret_cast<double*
>(mmap(
nullptr,
sizeof(
double)*ISOSPEC_G_FACT_TABLE_SIZE, PROT_READ | PROT_WRITE, MAP_ANONYMOUS | MAP_PRIVATE, -1, 0));
38 ret =
reinterpret_cast<double*
>(calloc(ISOSPEC_G_FACT_TABLE_SIZE,
sizeof(
double)));
40 std::atexit(release_g_lfact_table);
44 double* g_lfact_table = alloc_lfact_table();
47 double RationalApproximation(
double t)
51 double c[] = {2.515517, 0.802853, 0.010328};
52 double d[] = {1.432788, 0.189269, 0.001308};
53 return t - ((c[2]*t + c[1])*t + c[0]) /
54 (((d[2]*t + d[1])*t + d[0])*t + 1.0);
57 double NormalCDFInverse(
double p)
60 return -RationalApproximation( sqrt(-2.0*log(p)) );
62 return RationalApproximation( sqrt(-2.0*log(1-p)) );
65 double NormalCDFInverse(
double p,
double mean,
double stdev)
67 return mean + stdev * NormalCDFInverse(p);
70 double NormalCDF(
double x,
double mean,
double stdev)
72 x = (x-mean)/stdev * 0.7071067811865476;
75 double a1 = 0.254829592;
76 double a2 = -0.284496736;
77 double a3 = 1.421413741;
78 double a4 = -1.453152027;
79 double a5 = 1.061405429;
89 double t = 1.0/(1.0 + p*x);
90 double y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x);
92 return 0.5*(1.0 + sign*y);
95 double NormalPDF(
double x,
double mean,
double stdev)
97 double two_variance = stdev * stdev * 2.0;
98 double delta = x-mean;
99 return exp( -delta*delta / two_variance ) / sqrt( two_variance * pi );
102 const double sqrt_pi = 1.772453850905516027298167483341145182798;
104 double LowerIncompleteGamma2(
int a,
double x)
107 double exp_minus_x = exp(-x);
111 base = 1 - exp_minus_x;
117 base = sqrt_pi * erf(sqrt(x));
124 base = base * current_s - pow(x, current_s) * exp_minus_x;
131 double InverseLowerIncompleteGamma2(
int a,
double x)
134 double p = tgamma(a);
139 double v = LowerIncompleteGamma2(a, s);
144 }
while((p-l)*1000.0 > p);
149 std::random_device random_dev;
150 std::mt19937 random_gen(random_dev());
151 std::uniform_real_distribution<double> stdunif(0.0, 1.0);
153 size_t rdvariate_binom(
size_t tries,
double succ_prob, std::mt19937& rgen)
155 if (succ_prob >= 1.0)
157 return IsoSpec::boost_binomial_distribution_variate(tries, succ_prob, rgen);