IsoSpec  1.95
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 #include <cmath>
10 #include <cstdlib>
11 #include "isoMath.h"
12 #include "platform.h"
13 
14 namespace IsoSpec
15 {
16 
17 
18 void release_g_lfact_table()
19 {
20 #if ISOSPEC_GOT_MMAN
21  munmap(g_lfact_table, ISOSPEC_G_FACT_TABLE_SIZE*sizeof(double));
22 #else
23  free(g_lfact_table);
24 #endif
25 }
26 
27 double* alloc_lfact_table()
28 {
29  double* ret;
30 # if ISOSPEC_GOT_MMAN
31  ret = reinterpret_cast<double*>(mmap(nullptr, sizeof(double)*ISOSPEC_G_FACT_TABLE_SIZE, PROT_READ | PROT_WRITE, MAP_ANONYMOUS | MAP_PRIVATE, -1, 0));
32 #else
33  ret = reinterpret_cast<double*>(calloc(ISOSPEC_G_FACT_TABLE_SIZE, sizeof(double)));
34 #endif
35  std::atexit(release_g_lfact_table);
36  return ret;
37 }
38 
39 double* g_lfact_table = alloc_lfact_table();
40 
41 
42 double RationalApproximation(double t)
43 {
44  // Abramowitz and Stegun formula 26.2.23.
45  // The absolute value of the error should be less than 4.5 e-4.
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);
50 }
51 
52 double NormalCDFInverse(double p)
53 {
54 
55  if (p < 0.5)
56  return -RationalApproximation( sqrt(-2.0*log(p)) );
57  else
58  return RationalApproximation( sqrt(-2.0*log(1-p)) );
59 }
60 
61 double NormalCDFInverse(double p, double mean, double stdev)
62 {
63  return mean + stdev * NormalCDFInverse(p);
64 }
65 
66 double NormalCDF(double x, double mean, double stdev)
67 {
68  x = (x-mean)/stdev * 0.7071067811865476;
69 
70  // constants
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;
76  double p = 0.3275911;
77 
78  // Save the sign of x
79  int sign = 1;
80  if (x < 0)
81  sign = -1;
82  x = fabs(x);
83 
84  // A&S formula 7.1.26
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);
87 
88  return 0.5*(1.0 + sign*y);
89 }
90 
91 double NormalPDF(double x, double mean, double stdev)
92 {
93  double two_variance = stdev * stdev * 2.0;
94  double delta = x-mean;
95  return exp( -delta*delta / two_variance ) / sqrt( two_variance * pi );
96 }
97 
98 const double sqrt_pi = 1.772453850905516027298167483341145182798;
99 
100 double LowerIncompleteGamma2(int a, double x)
101 {
102  double base;
103  double exp_minus_x = exp(-x);
104  double current_s;
105  if(a % 2 == 0)
106  {
107  base = 1 - exp_minus_x;
108  current_s = 1.0;
109  a--;
110  }
111  else
112  {
113  base = sqrt_pi * erf(sqrt(x));
114  current_s = 0.5;
115  }
116 
117  a = a/2;
118  for(; a; a--)
119  {
120  base = base * current_s - pow(x, current_s) * exp_minus_x;
121  current_s += 1.0;
122  }
123 
124  return base;
125 }
126 
127 double InverseLowerIncompleteGamma2(int a, double x)
128 {
129  double l = 0.0;
130  double p = tgamma(a);
131  double s, v;
132 
133  do {
134  s = (l+p) / 2.0;
135  v = LowerIncompleteGamma2(a, s);
136  if (x < v)
137  p = s;
138  else
139  l = s;
140  } while((p-l)*1000.0>p);
141 
142  return s;
143 }
144 
145 
146 
147 } // namespace IsoSpec
148 
IsoSpec
Definition: allocator.cpp:21