| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 // $Id: Gamma.cc 68057 2013-03-13 14:46:00Z gcosmo $
26 // 27 //
27 // 28 //
28 // ------------------------------------------- 29 // ------------------------------------------------------------
29 // GEANT 4 class implementation 30 // GEANT 4 class implementation
30 // ------------------------------------------- 31 // ------------------------------------------------------------
31 32
32 #include <cmath> 33 #include <cmath>
33 #include <string.h> 34 #include <string.h>
34 #include "Gamma.hh" 35 #include "Gamma.hh"
35 36
36 MyGamma::MyGamma(){} 37 MyGamma::MyGamma(){}
37 38
38 MyGamma::~MyGamma(){} 39 MyGamma::~MyGamma(){}
39 40
40 //____________________________________________ 41 //____________________________________________________________________________
41 double MyGamma::Gamma(double z) 42 double MyGamma::Gamma(double z)
42 { 43 {
43 if (z <= 0) << 44 // Computation of gamma(z) for all z>0.
44 return 0; << 45 //
45 << 46 // The algorithm is based on the article by C.Lanczos [1] as denoted in
46 return std::tgamma(z); << 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);
47 } 57 }
48 58
49 //____________________________________________ 59 //____________________________________________________________________________
50 double MyGamma::Gamma(double a, double x) << 60 double MyGamma::Gamma(double a,double x)
51 { 61 {
52 // Computation of the incomplete gamma funct 62 // Computation of the incomplete gamma function P(a,x)
53 // 63 //
54 // The algorithm is based on the formulas an 64 // The algorithm is based on the formulas and code as denoted in
55 // Numerical Recipes 2nd ed. on p. 210-212 ( 65 // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
56 // 66 //
57 //--- Nve 14-nov-1998 UU-SAP Utrecht 67 //--- Nve 14-nov-1998 UU-SAP Utrecht
58 << 68
59 if (a <= 0 || x <= 0) return 0; 69 if (a <= 0 || x <= 0) return 0;
60 << 70
61 if (x < (a + 1)) << 71 if (x < (a+1)) return GamSer(a,x);
62 return GamSer(a, x); << 72 else return GamCf(a,x);
63 else <<
64 return GamCf(a, x); <<
65 } 73 }
66 74
67 //____________________________________________ 75 //____________________________________________________________________________
68 double MyGamma::GamCf(double a, double x) << 76 double MyGamma::GamCf(double a,double x)
69 { 77 {
70 // Computation of the incomplete gamma funct 78 // Computation of the incomplete gamma function P(a,x)
71 // via its continued fraction representation 79 // via its continued fraction representation.
72 // 80 //
73 // The algorithm is based on the formulas an 81 // The algorithm is based on the formulas and code as denoted in
74 // Numerical Recipes 2nd ed. on p. 210-212 ( 82 // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
75 // 83 //
76 //--- Nve 14-nov-1998 UU-SAP Utrecht 84 //--- Nve 14-nov-1998 UU-SAP Utrecht
77 << 85
78 int itmax = 100; // Maximum number of itera << 86 int itmax = 100; // Maximum number of iterations
79 double eps = 3.e-7; // Relative accuracy << 87 double eps = 3.e-7; // Relative accuracy
80 double fpmin = 1.e-30; // Smallest double v << 88 double fpmin = 1.e-30; // Smallest double value allowed here
81 << 89
82 if (a <= 0 || x <= 0) return 0; 90 if (a <= 0 || x <= 0) return 0;
83 << 91
84 double gln = LnGamma(a); 92 double gln = LnGamma(a);
85 double b = x + 1 - a; << 93 double b = x+1-a;
86 double c = 1 / fpmin; << 94 double c = 1/fpmin;
87 double d = 1 / b; << 95 double d = 1/b;
88 double h = d; << 96 double h = d;
89 double an, del; << 97 double an,del;
90 for (int i = 1; i <= itmax; i++) { << 98 for (int i=1; i<=itmax; i++) {
91 an = double(-i) * (double(i) - a); << 99 an = double(-i)*(double(i)-a);
92 b += 2; 100 b += 2;
93 d = an * d + b; << 101 d = an*d+b;
94 if (Abs(d) < fpmin) d = fpmin; 102 if (Abs(d) < fpmin) d = fpmin;
95 c = b + an / c; << 103 c = b+an/c;
96 if (Abs(c) < fpmin) c = fpmin; 104 if (Abs(c) < fpmin) c = fpmin;
97 d = 1 / d; << 105 d = 1/d;
98 del = d * c; << 106 del = d*c;
99 h = h * del; << 107 h = h*del;
100 if (Abs(del - 1) < eps) break; << 108 if (Abs(del-1) < eps) break;
101 // if (i==itmax) cout << "*GamCf(a,x)* a t << 109 //if (i==itmax) cout << "*GamCf(a,x)* a too large or itmax too small" << endl;
102 } 110 }
103 double v = Exp(-x + a * Log(x) - gln) * h; << 111 double v = Exp(-x+a*Log(x)-gln)*h;
104 return (1 - v); << 112 return (1-v);
105 } 113 }
106 114
107 //____________________________________________ 115 //____________________________________________________________________________
108 double MyGamma::GamSer(double a, double x) << 116 double MyGamma::GamSer(double a,double x)
109 { 117 {
110 // Computation of the incomplete gamma funct 118 // Computation of the incomplete gamma function P(a,x)
111 // via its series representation. 119 // via its series representation.
112 // 120 //
113 // The algorithm is based on the formulas an 121 // The algorithm is based on the formulas and code as denoted in
114 // Numerical Recipes 2nd ed. on p. 210-212 ( 122 // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
115 // 123 //
116 //--- Nve 14-nov-1998 UU-SAP Utrecht 124 //--- Nve 14-nov-1998 UU-SAP Utrecht
117 << 125
118 int itmax = 100; // Maximum number of itera << 126 int itmax = 100; // Maximum number of iterations
119 double eps = 3.e-7; // Relative accuracy << 127 double eps = 3.e-7; // Relative accuracy
120 << 128
121 if (a <= 0 || x <= 0) return 0; 129 if (a <= 0 || x <= 0) return 0;
122 << 130
123 double gln = LnGamma(a); 131 double gln = LnGamma(a);
124 double ap = a; << 132 double ap = a;
125 double sum = 1 / a; << 133 double sum = 1/a;
126 double del = sum; 134 double del = sum;
127 for (int n = 1; n <= itmax; n++) { << 135 for (int n=1; n<=itmax; n++) {
128 ap += 1; << 136 ap += 1;
129 del = del * x / ap; << 137 del = del*x/ap;
130 sum += del; 138 sum += del;
131 if (MyGamma::Abs(del) < Abs(sum * eps)) br << 139 if (MyGamma::Abs(del) < Abs(sum*eps)) break;
132 // if (n==itmax) cout << "*GamSer(a,x)* a << 140 //if (n==itmax) cout << "*GamSer(a,x)* a too large or itmax too small" << endl;
133 } 141 }
134 double v = sum * Exp(-x + a * Log(x) - gln); << 142 double v = sum*Exp(-x+a*Log(x)-gln);
135 return v; 143 return v;
136 } 144 }
137 145
>> 146
138 double MyGamma::LnGamma(double z) 147 double MyGamma::LnGamma(double z)
139 { 148 {
140 if (z <= 0) << 149 // Computation of ln[gamma(z)] for all z>0.
141 return 0; << 150 //
142 << 151 // The algorithm is based on the article by C.Lanczos [1] as denoted in
143 return std::lgamma(z); << 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;
144 } 178 }
145 179