IsoSpec  2.1.2
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 
19 namespace IsoSpec
20 {
21 
22 
23 void 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 
32 double* 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 
44 double* g_lfact_table = alloc_lfact_table();
45 
46 
47 double 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 
57 double 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 
65 double NormalCDFInverse(double p, double mean, double stdev)
66 {
67  return mean + stdev * NormalCDFInverse(p);
68 }
69 
70 double 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 
95 double 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 
102 const double sqrt_pi = 1.772453850905516027298167483341145182798;
103 
104 double 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 
131 double 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 
149 std::random_device random_dev;
150 std::mt19937 random_gen(random_dev());
151 std::uniform_real_distribution<double> stdunif(0.0, 1.0);
152 
153 size_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 
IsoSpec
Definition: allocator.cpp:20