VERB_code_2.3
erf.cpp
Go to the documentation of this file.
1 
8 #include <string>
9 #include "erf.h"
10 #include <iostream>
11 
12 using namespace std;
13 
14 bool str2bool(string str);
15 string bool2str(bool b);
16 
17 void nrerror(const char * msg) {
18  cout << msg;
19  return;
20 }
21 
23 double gammln(double xx)
24 {
25  //Internal arithmetic will be done in double precision, a nicety that you can omit if five-figure accuracy is good enough.
26  double x,y,tmp,ser;
27  static double cof[6]={76.18009172947146,-86.50532032941677, 24.01409824083091,-1.231739572450155,0.1208650973866179e-2,-0.5395239384953e-5};
28  int j;
29  y=x=xx;
30  tmp=x+5.5;
31  tmp -= (x+0.5)*log(tmp);
32  ser=1.000000000190015;
33  for (j=0;j<=5;j++) ser += cof[j]/++y;
34  return -tmp+log(2.5066282746310005*ser/x);
35 }
36 
37 
39 double gammp(double a, double x)
40 {
41  void gcf(double *gammcf, double a, double x, double *gln);
42  void gser(double *gamser, double a, double x, double *gln);
43 // void nrerror(char error_text[]);
44  double gamser,gammcf,gln;
45  if (x < 0.0 || a <= 0.0) nrerror("Invalid arguments in routine gammp");
46  if (x < (a+1.0)) { //Use the series representation.
47  gser(&gamser,a,x,&gln);
48  return gamser;
49  } else { //Use the continued fraction representation
50  gcf(&gammcf,a,x,&gln);
51  return 1.0-gammcf; //and take its complement.
52  }
53 }
54 
56 double gammq(double a, double x)
57 {
58  void gcf(double *gammcf, double a, double x, double *gln);
59  void gser(double *gamser, double a, double x, double *gln);
60 // void nrerror(char error_text[]);
61  double gamser,gammcf,gln;
62  if (x < 0.0 || a <= 0.0) nrerror("Invalid arguments in routine gammq");
63  if (x < (a+1.0)) { //Use the series representation
64  gser(&gamser,a,x,&gln);
65  return 1.0-gamser; //and take its complement.
66  } else { //Use the continued fraction representation.
67  gcf(&gammcf,a,x,&gln);
68  return gammcf;
69  }
70 }
71 
74 void gser(double *gamser, double a, double x, double *gln)
75 {
76  double gammln(double xx);
77 // void nrerror(char error_text[]);
78  int n;
79  double sum,del,ap;
80  *gln=gammln(a);
81  if (x <= 0.0) {
82  if (x < 0.0) nrerror("x less than 0 in routine gser");
83  *gamser=0.0;
84  return;
85  } else {
86  ap=a;
87  del=sum=1.0/a;
88  for (n=1;n<=ITMAX;n++) {
89  ++ap;
90  del *= x/ap;
91  sum += del;
92  if (fabs(del) < fabs(sum)*EPS) {
93  *gamser=sum*exp(-x+a*log(x)-(*gln));
94  return;
95  }
96  }
97  nrerror("a too large, ITMAX too small in routine gser");
98  return;
99  }
100 }
101 
103 void gcf(double *gammcf, double a, double x, double *gln)
104 {
105  double gammln(double xx);
106  double an;
107  double b;
108  double c1;
109  double d;
110  double del;
111  double h;
112 // void nrerror(char error_text[]);
113  int i;
114  *gln=gammln(a);
115  b=x+1.0-a; //Set up for evaluating continued fraction by modified Lentz's method (5.2) with b0 = 0.
116  c1=1.0/FPMIN;
117  d=1.0/b;
118  h=d;
119  for (i=1;i<=ITMAX;i++) { //Iterate to convergence.
120  an = -i*(i-a);
121  b += 2.0;
122  d=an*d+b;
123  if (fabs(d) < FPMIN) d=FPMIN;
124  c1=b+an/c1;
125  if (fabs(c1) < FPMIN) c1=FPMIN;
126  d=1.0/d;
127  del=d*c1;
128  h *= del;
129  if (fabs(del-1.0) < EPS) break;
130  }
131  if (i > ITMAX) nrerror("a too large, ITMAX too small in gcf");
132  *gammcf=exp(-x+a*log(x)-(*gln))*h; //Put factors in front.
133 }
134 
135 
137 double erf(double x)
138 {
139  double gammp(double a, double x);
140  return x < 0.0 ? -gammp(0.5,x*x) : gammp(0.5,x*x);
141 }
#define EPS
Relative accuracy.
Definition: erf.h:13
#define FPMIN
Number near the smallest representable doubleing-point number.
Definition: erf.h:15
double erf(double x)
Returns the error function erf(x).
Definition: erf.cpp:137
double gammln(double xx)
Returns the value ln[?(xx)] for xx > 0.
Definition: erf.cpp:23
static const double exp
exponential
bool str2bool(string str)
Definition: Parameters.cpp:583
double gammq(double a, double x)
Returns the incomplete gamma function Q(a, x) ? 1 ? P(a, x).
Definition: erf.cpp:56
#define ITMAX
Maximum allowed number of iterations.
Definition: erf.h:11
string bool2str(bool b)
Definition: Parameters.cpp:591
void gcf(double *gammcf, double a, double x, double *gln)
Returns the incomplete gamma function Q(a, x) evaluated by its continued fraction representation as g...
Definition: erf.cpp:103
void gser(double *gamser, double a, double x, double *gln)
Definition: erf.cpp:74
double gammp(double a, double x)
Returns the incomplete gamma function P(a, x).
Definition: erf.cpp:39