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 const double pi = 3.14159265358979323846264338328;
18 
19 void release_g_lfact_table()
20 {
21 #if ISOSPEC_GOT_MMAN
22  munmap(g_lfact_table, ISOSPEC_G_FACT_TABLE_SIZE*sizeof(double));
23 #else
24  free(g_lfact_table);
25 #endif
26 }
27 
28 double* alloc_lfact_table()
29 {
30  double* ret;
31 # if ISOSPEC_GOT_MMAN
32  ret = reinterpret_cast<double*>(mmap(nullptr, sizeof(double)*ISOSPEC_G_FACT_TABLE_SIZE, PROT_READ | PROT_WRITE, MAP_ANONYMOUS | MAP_PRIVATE, -1, 0));
33 #else
34  ret = reinterpret_cast<double*>(calloc(ISOSPEC_G_FACT_TABLE_SIZE, sizeof(double)));
35 #endif
36  std::atexit(release_g_lfact_table);
37  return ret;
38 }
39 
40 double* g_lfact_table = alloc_lfact_table();
41 
42 
43 double RationalApproximation(double t)
44 {
45  // Abramowitz and Stegun formula 26.2.23.
46  // The absolute value of the error should be less than 4.5 e-4.
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);
51 }
52 
53 double NormalCDFInverse(double p)
54 {
55 
56  if (p < 0.5)
57  return -RationalApproximation( sqrt(-2.0*log(p)) );
58  else
59  return RationalApproximation( sqrt(-2.0*log(1-p)) );
60 }
61 
62 double NormalCDFInverse(double p, double mean, double stdev)
63 {
64  return mean + stdev * NormalCDFInverse(p);
65 }
66 
67 double NormalCDF(double x, double mean, double stdev)
68 {
69  x = (x-mean)/stdev * 0.7071067811865476;
70 
71  // constants
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;
77  double p = 0.3275911;
78 
79  // Save the sign of x
80  int sign = 1;
81  if (x < 0)
82  sign = -1;
83  x = fabs(x);
84 
85  // A&S formula 7.1.26
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);
88 
89  return 0.5*(1.0 + sign*y);
90 }
91 
92 double NormalPDF(double x, double mean, double stdev)
93 {
94  double two_variance = stdev * stdev * 2.0;
95  double delta = x-mean;
96  return exp( -delta*delta / two_variance ) / sqrt( two_variance * pi );
97 }
98 
99 } // namespace IsoSpec
100 
IsoSpec
Definition: allocator.cpp:22