| 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,v 1.6 2006/06/29 19:14:28 gunter Exp $
>> 27 // GEANT4 tag $Name: geant4-08-01-patch-01 $
26 // 28 //
27 // 29 //
28 // ------------------------------------------- 30 // ------------------------------------------------------------
29 // GEANT 4 class implementation 31 // GEANT 4 class implementation
30 // ------------------------------------------- 32 // ------------------------------------------------------------
31 33
32 #include <cmath> 34 #include <cmath>
33 #include <string.h> 35 #include <string.h>
34 #include "Gamma.hh" 36 #include "Gamma.hh"
35 37
36 MyGamma::MyGamma(){} 38 MyGamma::MyGamma(){}
37 39
38 MyGamma::~MyGamma(){} 40 MyGamma::~MyGamma(){}
39 41
40 //____________________________________________ 42 //____________________________________________________________________________
41 double MyGamma::Gamma(double z) 43 double MyGamma::Gamma(double z)
42 { 44 {
43 if (z <= 0) << 45 // Computation of gamma(z) for all z>0.
44 return 0; << 46 //
45 << 47 // The algorithm is based on the article by C.Lanczos [1] as denoted in
46 return std::tgamma(z); << 48 // Numerical Recipes 2nd ed. on p. 207 (W.H.Press et al.).
>> 49 //
>> 50 // [1] C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86.
>> 51 //
>> 52 //--- Nve 14-nov-1998 UU-SAP Utrecht
>> 53
>> 54 if (z<=0) return 0;
>> 55
>> 56 double v = LnGamma(z);
>> 57 return std::exp(v);
47 } 58 }
48 59
49 //____________________________________________ 60 //____________________________________________________________________________
50 double MyGamma::Gamma(double a, double x) << 61 double MyGamma::Gamma(double a,double x)
51 { 62 {
52 // Computation of the incomplete gamma funct 63 // Computation of the incomplete gamma function P(a,x)
53 // 64 //
54 // The algorithm is based on the formulas an 65 // The algorithm is based on the formulas and code as denoted in
55 // Numerical Recipes 2nd ed. on p. 210-212 ( 66 // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
56 // 67 //
57 //--- Nve 14-nov-1998 UU-SAP Utrecht 68 //--- Nve 14-nov-1998 UU-SAP Utrecht
58 << 69
59 if (a <= 0 || x <= 0) return 0; 70 if (a <= 0 || x <= 0) return 0;
60 << 71
61 if (x < (a + 1)) << 72 if (x < (a+1)) return GamSer(a,x);
62 return GamSer(a, x); << 73 else return GamCf(a,x);
63 else <<
64 return GamCf(a, x); <<
65 } 74 }
66 75
67 //____________________________________________ 76 //____________________________________________________________________________
68 double MyGamma::GamCf(double a, double x) << 77 double MyGamma::GamCf(double a,double x)
69 { 78 {
70 // Computation of the incomplete gamma funct 79 // Computation of the incomplete gamma function P(a,x)
71 // via its continued fraction representation 80 // via its continued fraction representation.
72 // 81 //
73 // The algorithm is based on the formulas an 82 // The algorithm is based on the formulas and code as denoted in
74 // Numerical Recipes 2nd ed. on p. 210-212 ( 83 // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
75 // 84 //
76 //--- Nve 14-nov-1998 UU-SAP Utrecht 85 //--- Nve 14-nov-1998 UU-SAP Utrecht
77 << 86
78 int itmax = 100; // Maximum number of itera << 87 int itmax = 100; // Maximum number of iterations
79 double eps = 3.e-7; // Relative accuracy << 88 double eps = 3.e-7; // Relative accuracy
80 double fpmin = 1.e-30; // Smallest double v << 89 double fpmin = 1.e-30; // Smallest double value allowed here
81 << 90
82 if (a <= 0 || x <= 0) return 0; 91 if (a <= 0 || x <= 0) return 0;
83 << 92
84 double gln = LnGamma(a); 93 double gln = LnGamma(a);
85 double b = x + 1 - a; << 94 double b = x+1-a;
86 double c = 1 / fpmin; << 95 double c = 1/fpmin;
87 double d = 1 / b; << 96 double d = 1/b;
88 double h = d; << 97 double h = d;
89 double an, del; << 98 double an,del;
90 for (int i = 1; i <= itmax; i++) { << 99 for (int i=1; i<=itmax; i++) {
91 an = double(-i) * (double(i) - a); << 100 an = double(-i)*(double(i)-a);
92 b += 2; 101 b += 2;
93 d = an * d + b; << 102 d = an*d+b;
94 if (Abs(d) < fpmin) d = fpmin; 103 if (Abs(d) < fpmin) d = fpmin;
95 c = b + an / c; << 104 c = b+an/c;
96 if (Abs(c) < fpmin) c = fpmin; 105 if (Abs(c) < fpmin) c = fpmin;
97 d = 1 / d; << 106 d = 1/d;
98 del = d * c; << 107 del = d*c;
99 h = h * del; << 108 h = h*del;
100 if (Abs(del - 1) < eps) break; << 109 if (Abs(del-1) < eps) break;
101 // if (i==itmax) cout << "*GamCf(a,x)* a t << 110 //if (i==itmax) cout << "*GamCf(a,x)* a too large or itmax too small" << endl;
102 } 111 }
103 double v = Exp(-x + a * Log(x) - gln) * h; << 112 double v = Exp(-x+a*Log(x)-gln)*h;
104 return (1 - v); << 113 return (1-v);
105 } 114 }
106 115
107 //____________________________________________ 116 //____________________________________________________________________________
108 double MyGamma::GamSer(double a, double x) << 117 double MyGamma::GamSer(double a,double x)
109 { 118 {
110 // Computation of the incomplete gamma funct 119 // Computation of the incomplete gamma function P(a,x)
111 // via its series representation. 120 // via its series representation.
112 // 121 //
113 // The algorithm is based on the formulas an 122 // The algorithm is based on the formulas and code as denoted in
114 // Numerical Recipes 2nd ed. on p. 210-212 ( 123 // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
115 // 124 //
116 //--- Nve 14-nov-1998 UU-SAP Utrecht 125 //--- Nve 14-nov-1998 UU-SAP Utrecht
117 << 126
118 int itmax = 100; // Maximum number of itera << 127 int itmax = 100; // Maximum number of iterations
119 double eps = 3.e-7; // Relative accuracy << 128 double eps = 3.e-7; // Relative accuracy
120 << 129
121 if (a <= 0 || x <= 0) return 0; 130 if (a <= 0 || x <= 0) return 0;
122 << 131
123 double gln = LnGamma(a); 132 double gln = LnGamma(a);
124 double ap = a; << 133 double ap = a;
125 double sum = 1 / a; << 134 double sum = 1/a;
126 double del = sum; 135 double del = sum;
127 for (int n = 1; n <= itmax; n++) { << 136 for (int n=1; n<=itmax; n++) {
128 ap += 1; << 137 ap += 1;
129 del = del * x / ap; << 138 del = del*x/ap;
130 sum += del; 139 sum += del;
131 if (MyGamma::Abs(del) < Abs(sum * eps)) br << 140 if (MyGamma::Abs(del) < Abs(sum*eps)) break;
132 // if (n==itmax) cout << "*GamSer(a,x)* a << 141 //if (n==itmax) cout << "*GamSer(a,x)* a too large or itmax too small" << endl;
133 } 142 }
134 double v = sum * Exp(-x + a * Log(x) - gln); << 143 double v = sum*Exp(-x+a*Log(x)-gln);
135 return v; 144 return v;
136 } 145 }
137 146
>> 147
138 double MyGamma::LnGamma(double z) 148 double MyGamma::LnGamma(double z)
139 { 149 {
140 if (z <= 0) << 150 // Computation of ln[gamma(z)] for all z>0.
141 return 0; << 151 //
142 << 152 // The algorithm is based on the article by C.Lanczos [1] as denoted in
143 return std::lgamma(z); << 153 // Numerical Recipes 2nd ed. on p. 207 (W.H.Press et al.).
>> 154 //
>> 155 // [1] C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86.
>> 156 //
>> 157 // The accuracy of the result is better than 2e-10.
>> 158 //
>> 159 //--- Nve 14-nov-1998 UU-SAP Utrecht
>> 160
>> 161 if (z<=0) return 0;
>> 162
>> 163 // Coefficients for the series expansion
>> 164 double c[7] = { 2.5066282746310005, 76.18009172947146, -86.50532032941677
>> 165 ,24.01409824083091, -1.231739572450155, 0.1208650973866179e-2
>> 166 ,-0.5395239384953e-5};
>> 167
>> 168 double x = z;
>> 169 double y = x;
>> 170 double tmp = x+5.5;
>> 171 tmp = (x+0.5)*Log(tmp)-tmp;
>> 172 double ser = 1.000000000190015;
>> 173 for (int i=1; i<7; i++) {
>> 174 y += 1;
>> 175 ser += c[i]/y;
>> 176 }
>> 177 double v = tmp+Log(c[0]*ser/x);
>> 178 return v;
144 } 179 }
145 180