Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Gamma.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: Gamma.cc 69796 2013-05-15 13:26:12Z gcosmo $
27 //
28 //
29 // ------------------------------------------------------------
30 // GEANT 4 class implementation
31 // ------------------------------------------------------------
32 
33 #include <cmath>
34 #include <string.h>
35 #include "Gamma.hh"
36 
38 
40 
41 //____________________________________________________________________________
42 double MyGamma::Gamma(double z)
43 {
44  // Computation of gamma(z) for all z>0.
45  //
46  // The algorithm is based on the article by C.Lanczos [1] as denoted in
47  // Numerical Recipes 2nd ed. on p. 207 (W.H.Press et al.).
48  //
49  // [1] C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86.
50  //
51  //--- Nve 14-nov-1998 UU-SAP Utrecht
52 
53  if (z<=0) return 0;
54 
55  double v = LnGamma(z);
56  return std::exp(v);
57 }
58 
59 //____________________________________________________________________________
60 double MyGamma::Gamma(double a,double x)
61 {
62  // Computation of the incomplete gamma function P(a,x)
63  //
64  // The algorithm is based on the formulas and code as denoted in
65  // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
66  //
67  //--- Nve 14-nov-1998 UU-SAP Utrecht
68 
69  if (a <= 0 || x <= 0) return 0;
70 
71  if (x < (a+1)) return GamSer(a,x);
72  else return GamCf(a,x);
73 }
74 
75 //____________________________________________________________________________
76 double MyGamma::GamCf(double a,double x)
77 {
78  // Computation of the incomplete gamma function P(a,x)
79  // via its continued fraction representation.
80  //
81  // The algorithm is based on the formulas and code as denoted in
82  // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
83  //
84  //--- Nve 14-nov-1998 UU-SAP Utrecht
85 
86  int itmax = 100; // Maximum number of iterations
87  double eps = 3.e-7; // Relative accuracy
88  double fpmin = 1.e-30; // Smallest double value allowed here
89 
90  if (a <= 0 || x <= 0) return 0;
91 
92  double gln = LnGamma(a);
93  double b = x+1-a;
94  double c = 1/fpmin;
95  double d = 1/b;
96  double h = d;
97  double an,del;
98  for (int i=1; i<=itmax; i++) {
99  an = double(-i)*(double(i)-a);
100  b += 2;
101  d = an*d+b;
102  if (Abs(d) < fpmin) d = fpmin;
103  c = b+an/c;
104  if (Abs(c) < fpmin) c = fpmin;
105  d = 1/d;
106  del = d*c;
107  h = h*del;
108  if (Abs(del-1) < eps) break;
109  //if (i==itmax) cout << "*GamCf(a,x)* a too large or itmax too small" << endl;
110  }
111  double v = Exp(-x+a*Log(x)-gln)*h;
112  return (1-v);
113 }
114 
115 //____________________________________________________________________________
116 double MyGamma::GamSer(double a,double x)
117 {
118  // Computation of the incomplete gamma function P(a,x)
119  // via its series representation.
120  //
121  // The algorithm is based on the formulas and code as denoted in
122  // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
123  //
124  //--- Nve 14-nov-1998 UU-SAP Utrecht
125 
126  int itmax = 100; // Maximum number of iterations
127  double eps = 3.e-7; // Relative accuracy
128 
129  if (a <= 0 || x <= 0) return 0;
130 
131  double gln = LnGamma(a);
132  double ap = a;
133  double sum = 1/a;
134  double del = sum;
135  for (int n=1; n<=itmax; n++) {
136  ap += 1;
137  del = del*x/ap;
138  sum += del;
139  if (MyGamma::Abs(del) < Abs(sum*eps)) break;
140  //if (n==itmax) cout << "*GamSer(a,x)* a too large or itmax too small" << endl;
141  }
142  double v = sum*Exp(-x+a*Log(x)-gln);
143  return v;
144 }
145 
146 
147 double MyGamma::LnGamma(double z)
148 {
149  // Computation of ln[gamma(z)] for all z>0.
150  //
151  // The algorithm is based on the article by C.Lanczos [1] as denoted in
152  // Numerical Recipes 2nd ed. on p. 207 (W.H.Press et al.).
153  //
154  // [1] C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86.
155  //
156  // The accuracy of the result is better than 2e-10.
157  //
158  //--- Nve 14-nov-1998 UU-SAP Utrecht
159 
160  if (z<=0) return 0;
161 
162  // Coefficients for the series expansion
163  double c[7] = { 2.5066282746310005, 76.18009172947146, -86.50532032941677
164  ,24.01409824083091, -1.231739572450155, 0.1208650973866179e-2
165  ,-0.5395239384953e-5};
166 
167  double x = z;
168  double y = x;
169  double tmp = x+5.5;
170  tmp = (x+0.5)*Log(tmp)-tmp;
171  double ser = 1.000000000190015;
172  for (int i=1; i<7; i++) {
173  y += 1;
174  ser += c[i]/y;
175  }
176  double v = tmp+Log(c[0]*ser/x);
177  return v;
178 }