LCFIVertex  0.7.2
util.cpp
1 #include "../inc/util.h"
2 #include <math.h>
3 
4 namespace vertex_lcfi { namespace util{
5 
6 double gamma( double a, double x)
7  {
8  if (a <= 0 || x <= 0) return 0;
9 
10  if (x < (a+1)) return gamSer(a,x);
11  else return gamCf(a,x);
12  }
13 
14  double gamCf(double a,double x)
15  {
16  int itmax = 100; // Maximum number of iterations
17  double eps = 3.e-14; // Relative accuracy
18  double fpmin = 1.e-30; // Smallest double value allowed here
19 
20  if (a <= 0 || x <= 0) return 0;
21 
22  double gln = lnGamma(a);
23  double b = x+1-a;
24  double c = 1/fpmin;
25  double d = 1/b;
26  double h = d;
27  double an,del;
28  for (int i=1; i<=itmax; i++) {
29  an = double(-i)*(double(i)-a);
30  b += 2;
31  d = an*d+b;
32  if (fabs(d) < fpmin) d = fpmin;
33  c = b+an/c;
34  if (fabs(c) < fpmin) c = fpmin;
35  d = 1/d;
36  del = d*c;
37  h = h*del;
38  if (fabs(del-1) < eps) break;
39  //if (i==itmax) cout << "*GamCf(a,x)* a too large or itmax too small" << endl;
40  }
41  double v = exp(-x+a*log(x)-gln)*h;
42  return (1-v);
43  }
44 
45 
46  double gamSer(double a,double x)
47  {
48  int itmax = 100; // Maximum number of iterations
49  double eps = 3.e-14; // Relative accuracy
50 
51  if (a <= 0 || x <= 0) return 0;
52 
53  double gln = lnGamma(a);
54  double ap = a;
55  double sum = 1/a;
56  double del = sum;
57  for (int n=1; n<=itmax; n++) {
58  ap += 1;
59  del = del*x/ap;
60  sum += del;
61  if ( fabs(del) < fabs(sum*eps)) break;
62  //if (n==itmax) cout << "*GamSer(a,x)* a too large or itmax too small" << endl;
63  }
64  double v = sum*exp(-x+a*log(x)-gln);
65  return v;
66  }
67 
68  double lnGamma(double z)
69  {
70  if (z<=0) return 0;
71 
72  // Coefficients for the series expansion
73  double c[7] = { 2.5066282746310005, 76.18009172947146, -86.50532032941677
74  ,24.01409824083091, -1.231739572450155, 0.1208650973866179e-2
75  ,-0.5395239384953e-5};
76 
77  double x = z;
78  double y = x;
79  double tmp = x+5.5;
80  tmp = (x+0.5)*log(tmp)-tmp;
81  double ser = 1.000000000190015;
82  for (int i=1; i<7; i++) {
83  y += 1;
84  ser += c[i]/y;
85  }
86  double v = tmp+log(c[0]*ser/x);
87  return v;
88  }
89 
90 }}
double gamCf(double a, double x)
Computation of the incomplete gamma function P(a,x) via its continued fraction representation.
Definition: util.cpp:14
double gamSer(double a, double x)
Incomplete gamma function P(a,x) via its series representation.
Definition: util.cpp:46
double lnGamma(double z)
Computation of ln[gamma(z)] for all z>0.
Definition: util.cpp:68
double gamma(double a, double x)
Computation of the upper incomplete gamma function P(a,x)
Definition: util.cpp:6