IsoSpec 2.2.1
isoMath.cpp
1/*
2 * This file has been released into public domain by John D. Cook
3 * and is used here with some slight modifications (which are hereby
4 * also released into public domain),
5 *
6 * This file is part of IsoSpec.
7 */
8
9
10// NOLINT(legal/copyright)
11
12
13#include <cmath>
14#include <cstdlib>
15#include "isoMath.h"
16#include "platform.h"
17#include "btrd.h"
18
19namespace IsoSpec
20{
21
22
23void release_g_lfact_table()
24{
25#if ISOSPEC_GOT_MMAN
26 munmap(g_lfact_table, ISOSPEC_G_FACT_TABLE_SIZE*sizeof(double));
27#else
28 free(g_lfact_table);
29#endif
30}
31
32double* alloc_lfact_table()
33{
34 double* ret;
35# if ISOSPEC_GOT_MMAN
36 ret = reinterpret_cast<double*>(mmap(nullptr, sizeof(double)*ISOSPEC_G_FACT_TABLE_SIZE, PROT_READ | PROT_WRITE, MAP_ANONYMOUS | MAP_PRIVATE, -1, 0));
37#else
38 ret = reinterpret_cast<double*>(calloc(ISOSPEC_G_FACT_TABLE_SIZE, sizeof(double)));
39#endif
40 std::atexit(release_g_lfact_table);
41 return ret;
42}
43
44double* g_lfact_table = alloc_lfact_table();
45
46
47double RationalApproximation(double t)
48{
49 // Abramowitz and Stegun formula 26.2.23.
50 // The absolute value of the error should be less than 4.5 e-4.
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);
55}
56
57double NormalCDFInverse(double p)
58{
59 if (p < 0.5)
60 return -RationalApproximation( sqrt(-2.0*log(p)) );
61 else
62 return RationalApproximation( sqrt(-2.0*log(1-p)) );
63}
64
65double NormalCDFInverse(double p, double mean, double stdev)
66{
67 return mean + stdev * NormalCDFInverse(p);
68}
69
70double NormalCDF(double x, double mean, double stdev)
71{
72 x = (x-mean)/stdev * 0.7071067811865476;
73
74 // constants
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;
80 double p = 0.3275911;
81
82 // Save the sign of x
83 int sign = 1;
84 if (x < 0)
85 sign = -1;
86 x = fabs(x);
87
88 // A&S formula 7.1.26
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);
91
92 return 0.5*(1.0 + sign*y);
93}
94
95double NormalPDF(double x, double mean, double stdev)
96{
97 double two_variance = stdev * stdev * 2.0;
98 double delta = x-mean;
99 return exp( -delta*delta / two_variance ) / sqrt( two_variance * pi );
100}
101
102const double sqrt_pi = 1.772453850905516027298167483341145182798;
103
104double LowerIncompleteGamma2(int a, double x)
105{
106 double base;
107 double exp_minus_x = exp(-x);
108 double current_s;
109 if(a % 2 == 0)
110 {
111 base = 1 - exp_minus_x;
112 current_s = 1.0;
113 a--;
114 }
115 else
116 {
117 base = sqrt_pi * erf(sqrt(x));
118 current_s = 0.5;
119 }
120
121 a = a/2;
122 for(; a; a--)
123 {
124 base = base * current_s - pow(x, current_s) * exp_minus_x;
125 current_s += 1.0;
126 }
127
128 return base;
129}
130
131double InverseLowerIncompleteGamma2(int a, double x)
132{
133 double l = 0.0;
134 double p = tgamma(a);
135 double s;
136
137 do {
138 s = (l+p) / 2.0;
139 double v = LowerIncompleteGamma2(a, s);
140 if (x < v)
141 p = s;
142 else
143 l = s;
144 } while((p-l)*1000.0 > p);
145
146 return s;
147}
148
149std::random_device random_dev;
150std::mt19937 random_gen(random_dev());
151std::uniform_real_distribution<double> stdunif(0.0, 1.0);
152
153size_t rdvariate_binom(size_t tries, double succ_prob, std::mt19937& rgen)
154{
155 if (succ_prob >= 1.0)
156 return tries;
157 return IsoSpec::boost_binomial_distribution_variate(tries, succ_prob, rgen);
158}
159
160
161
162} // namespace IsoSpec
163