| Geant4 Cross Reference |
1 /* gamma.c 1 /* gamma.c
2 * 2 *
3 * Gamma function 3 * Gamma function
4 * 4 *
5 * DESCRIPTION: 5 * DESCRIPTION:
6 * 6 *
7 * Returns gamma function of the argument. Th 7 * Returns gamma function of the argument. The result is
8 * correctly signed, and the sign (+1 or -1) i 8 * correctly signed, and the sign (+1 or -1) is also
9 * returned in a global (extern) variable name 9 * returned in a global (extern) variable named sgngam.
10 * This variable is also filled in by the loga 10 * This variable is also filled in by the logarithmic gamma
11 * function lgam(). 11 * function lgam().
12 * 12 *
13 * Arguments |x| <= 34 are reduced by recurren 13 * Arguments |x| <= 34 are reduced by recurrence and the function
14 * approximated by a rational function of degr 14 * approximated by a rational function of degree 6/7 in the
15 * interval (2,3). Large arguments are handle 15 * interval (2,3). Large arguments are handled by Stirling's
16 * formula. Large negative arguments are made 16 * formula. Large negative arguments are made positive using
17 * a reflection formula. 17 * a reflection formula.
18 * 18 *
19 * 19 *
20 * ACCURACY: 20 * ACCURACY:
21 * 21 *
22 * Relative error: 22 * Relative error:
23 * arithmetic domain # trials peak 23 * arithmetic domain # trials peak rms
24 * IEEE -170,-33 20000 2.3e-1 24 * IEEE -170,-33 20000 2.3e-15 3.3e-16
25 * IEEE -33, 33 20000 9.4e-1 25 * IEEE -33, 33 20000 9.4e-16 2.2e-16
26 * IEEE 33, 171.6 20000 2.3e-1 26 * IEEE 33, 171.6 20000 2.3e-15 3.2e-16
27 * 27 *
28 * Error for arguments outside the test range 28 * Error for arguments outside the test range will be larger
29 * owing to error amplification by the exponen 29 * owing to error amplification by the exponential function.
30 * 30 *
31 */ 31 */
32 /* lgam() 32 /* lgam()
33 * 33 *
34 * Natural logarithm of gamma function 34 * Natural logarithm of gamma function
35 * 35 *
36 * 36 *
37 * DESCRIPTION: 37 * DESCRIPTION:
38 * 38 *
39 * Returns the base e (2.718...) logarithm of 39 * Returns the base e (2.718...) logarithm of the absolute
40 * value of the gamma function of the argument 40 * value of the gamma function of the argument.
41 * The sign (+1 or -1) of the gamma function i 41 * The sign (+1 or -1) of the gamma function is returned in a
42 * global (extern) variable named sgngam. 42 * global (extern) variable named sgngam.
43 * 43 *
44 * For arguments greater than 13, the logarith 44 * For arguments greater than 13, the logarithm of the gamma
45 * function is approximated by the logarithmic 45 * function is approximated by the logarithmic version of
46 * Stirling's formula using a polynomial appro 46 * Stirling's formula using a polynomial approximation of
47 * degree 4. Arguments between -33 and +33 are 47 * degree 4. Arguments between -33 and +33 are reduced by
48 * recurrence to the interval [2,3] of a ratio 48 * recurrence to the interval [2,3] of a rational approximation.
49 * The cosecant reflection formula is employed 49 * The cosecant reflection formula is employed for arguments
50 * less than -33. 50 * less than -33.
51 * 51 *
52 * Arguments greater than MAXLGM return DBL_MA 52 * Arguments greater than MAXLGM return DBL_MAX and an error
53 * message. MAXLGM = 2.556348e305 for IEEE ar 53 * message. MAXLGM = 2.556348e305 for IEEE arithmetic.
54 * 54 *
55 * 55 *
56 * ACCURACY: 56 * ACCURACY:
57 * 57 *
58 * arithmetic domain # trials 58 * arithmetic domain # trials peak rms
59 * IEEE 0, 3 28000 5 59 * IEEE 0, 3 28000 5.4e-16 1.1e-16
60 * IEEE 2.718, 2.556e305 40000 3 60 * IEEE 2.718, 2.556e305 40000 3.5e-16 8.3e-17
61 * The error criterion was relative when the f 61 * The error criterion was relative when the function magnitude
62 * was greater than one but absolute when it w 62 * was greater than one but absolute when it was less than one.
63 * 63 *
64 * The following test used the relative error 64 * The following test used the relative error criterion, though
65 * at certain points the relative error could 65 * at certain points the relative error could be much higher than
66 * indicated. 66 * indicated.
67 * IEEE -200, -4 10000 4 67 * IEEE -200, -4 10000 4.8e-16 1.3e-16
68 * 68 *
69 */ 69 */
70 /* gamma.c */ 70 /* gamma.c */
71 /* gamma function */ 71 /* gamma function */
72 72
73 /* 73 /*
74 Cephes Math Library Release 2.8: June, 2000 74 Cephes Math Library Release 2.8: June, 2000
75 Copyright 1984, 1987, 1989, 1992, 2000 by Step 75 Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
76 */ 76 */
77 77
78 #include "nf_specialFunctions.h" 78 #include "nf_specialFunctions.h"
79 79
80 #if defined __cplusplus 80 #if defined __cplusplus
81 #include "G4Exp.hh" 81 #include "G4Exp.hh"
82 #include "G4Log.hh" 82 #include "G4Log.hh"
83 #include "G4Pow.hh" 83 #include "G4Pow.hh"
84 namespace GIDI { 84 namespace GIDI {
85 using namespace GIDI; 85 using namespace GIDI;
86 using namespace std; 86 using namespace std;
87 #endif 87 #endif
88 88
89 static double P[] = { 1.60119522476751861407E- 89 static double P[] = { 1.60119522476751861407E-4, 1.19135147006586384913E-3, 1.04213797561761569935E-2, 4.76367800457137231464E-2,
90 2.07448227648435975150E- 90 2.07448227648435975150E-1, 4.94214826801497100753E-1, 9.99999999999999996796E-1 };
91 static double Q[] = { -2.31581873324120129819E 91 static double Q[] = { -2.31581873324120129819E-5, 5.39605580493303397842E-4, -4.45641913851797240494E-3, 1.18139785222060435552E-2,
92 3.58236398605498653373E 92 3.58236398605498653373E-2, -2.34591795718243348568E-1, 7.14304917030273074085E-2, 1.00000000000000000320E0 };
93 #define MAXGAM 171.624376956302725 93 #define MAXGAM 171.624376956302725
94 static double LOGPI = 1.14472988584940017414; 94 static double LOGPI = 1.14472988584940017414;
95 static double SQTPI = 2.50662827463100050242E0 95 static double SQTPI = 2.50662827463100050242E0;
96 96
97 /* Stirling's formula for the gamma function * 97 /* Stirling's formula for the gamma function */
98 static double STIR[5] = { 7.873113957930936284 98 static double STIR[5] = { 7.873113957930936284e-4, -2.2954996161337812638e-4, -2.6813261780578123283e-3, 3.472222216054586673e-3, 8.3333333333348225713e-2 };
99 #define MAXSTIR 143.01608 99 #define MAXSTIR 143.01608
100 100
101 static double stirf( double x, nfu_status *sta 101 static double stirf( double x, nfu_status *status );
102 static double lgam( double x, int *sgngam, nfu 102 static double lgam( double x, int *sgngam, nfu_status *status );
103 /* 103 /*
104 ********************************************** 104 ************************************************************
105 */ 105 */
106 static double stirf( double x, nfu_status * /* 106 static double stirf( double x, nfu_status * /*status*/ ) {
107 /* Gamma function computed by Stirling's formu 107 /* Gamma function computed by Stirling's formula. The polynomial STIR is valid for 33 <= x <= 172. */
108 108
109 double y, w, v; 109 double y, w, v;
110 110
111 w = 1.0 / x; 111 w = 1.0 / x;
112 w = 1.0 + w * nf_polevl( w, STIR, 4 ); 112 w = 1.0 + w * nf_polevl( w, STIR, 4 );
113 y = G4Exp( x ); 113 y = G4Exp( x );
114 if( x > MAXSTIR ) { /* 114 if( x > MAXSTIR ) { /* Avoid overflow in pow() */
115 v = G4Pow::GetInstance()->powA( x, 0.5 115 v = G4Pow::GetInstance()->powA( x, 0.5 * x - 0.25 );
116 y = v * (v / y); } 116 y = v * (v / y); }
117 else { 117 else {
118 y = G4Pow::GetInstance()->powA( x, x - 118 y = G4Pow::GetInstance()->powA( x, x - 0.5 ) / y;
119 } 119 }
120 y = SQTPI * y * w; 120 y = SQTPI * y * w;
121 return( y ); 121 return( y );
122 } 122 }
123 /* 123 /*
124 ********************************************** 124 ************************************************************
125 */ 125 */
126 double nf_gammaFunction( double x, nfu_status 126 double nf_gammaFunction( double x, nfu_status *status ) {
127 127
128 double p, q, z; 128 double p, q, z;
129 int i, sgngam = 1; 129 int i, sgngam = 1;
130 130
131 *status = nfu_badInput; 131 *status = nfu_badInput;
132 if( !isfinite( x ) ) return( x ); 132 if( !isfinite( x ) ) return( x );
133 *status = nfu_Okay; 133 *status = nfu_Okay;
134 134
135 q = fabs( x ); 135 q = fabs( x );
136 136
137 if( q > 33.0 ) { 137 if( q > 33.0 ) {
138 if( x < 0.0 ) { 138 if( x < 0.0 ) {
139 p = floor( q ); 139 p = floor( q );
140 if( p == q ) goto goverf; 140 if( p == q ) goto goverf;
141 i = (int) p; 141 i = (int) p;
142 if( ( i & 1 ) == 0 ) sgngam = -1; 142 if( ( i & 1 ) == 0 ) sgngam = -1;
143 z = q - p; 143 z = q - p;
144 if( z > 0.5 ) { 144 if( z > 0.5 ) {
145 p += 1.0; 145 p += 1.0;
146 z = q - p; 146 z = q - p;
147 } 147 }
148 z = q * sin( M_PI * z ); 148 z = q * sin( M_PI * z );
149 if( z == 0.0 ) goto goverf; 149 if( z == 0.0 ) goto goverf;
150 z = M_PI / ( fabs( z ) * stirf( q, 150 z = M_PI / ( fabs( z ) * stirf( q, status ) );
151 } 151 }
152 else { 152 else {
153 z = stirf( x, status ); 153 z = stirf( x, status );
154 } 154 }
155 return( sgngam * z ); 155 return( sgngam * z );
156 } 156 }
157 157
158 z = 1.0; 158 z = 1.0;
159 while( x >= 3.0 ) { 159 while( x >= 3.0 ) {
160 x -= 1.0; 160 x -= 1.0;
161 z *= x; 161 z *= x;
162 } // Loop checking, 11.06.2015, T. Koi 162 } // Loop checking, 11.06.2015, T. Koi
163 163
164 while( x < 0.0 ) { 164 while( x < 0.0 ) {
165 if( x > -1.E-9 ) goto small; 165 if( x > -1.E-9 ) goto small;
166 z /= x; 166 z /= x;
167 x += 1.0; 167 x += 1.0;
168 } // Loop checking, 11.06.2015, T. Koi 168 } // Loop checking, 11.06.2015, T. Koi
169 169
170 while( x < 2.0 ) { 170 while( x < 2.0 ) {
171 if( x < 1.e-9 ) goto small; 171 if( x < 1.e-9 ) goto small;
172 z /= x; 172 z /= x;
173 x += 1.0; 173 x += 1.0;
174 } // Loop checking, 11.06.2015, T. Koi 174 } // Loop checking, 11.06.2015, T. Koi
175 175
176 if( x == 2.0 ) return( z ); 176 if( x == 2.0 ) return( z );
177 177
178 x -= 2.0; 178 x -= 2.0;
179 p = nf_polevl( x, P, 6 ); 179 p = nf_polevl( x, P, 6 );
180 q = nf_polevl( x, Q, 7 ); 180 q = nf_polevl( x, Q, 7 );
181 return( z * p / q ); 181 return( z * p / q );
182 182
183 small: 183 small:
184 if( x == 0.0 ) goto goverf; 184 if( x == 0.0 ) goto goverf;
185 return( z / ( ( 1.0 + 0.5772156649015329 * 185 return( z / ( ( 1.0 + 0.5772156649015329 * x ) * x ) );
186 186
187 goverf: 187 goverf:
188 return( sgngam * DBL_MAX ); 188 return( sgngam * DBL_MAX );
189 } 189 }
190 190
191 /* A[]: Stirling's formula expansion of log ga 191 /* A[]: Stirling's formula expansion of log gamma
192 * B[], C[]: log gamma function between 2 and 192 * B[], C[]: log gamma function between 2 and 3
193 */ 193 */
194 static double A[] = { 8.11614167470508450300E- 194 static double A[] = { 8.11614167470508450300E-4, -5.95061904284301438324E-4, 7.93650340457716943945E-4,
195 -2.77777777730099687205E- 195 -2.77777777730099687205E-3, 8.33333333333331927722E-2 };
196 static double B[] = { -1.37825152569120859100E 196 static double B[] = { -1.37825152569120859100E3, -3.88016315134637840924E4, -3.31612992738871184744E5,
197 -1.16237097492762307383E 197 -1.16237097492762307383E6, -1.72173700820839662146E6, -8.53555664245765465627E5 };
198 static double C[] = { -3.51815701436523470549E 198 static double C[] = { -3.51815701436523470549E2, -1.70642106651881159223E4, -2.20528590553854454839E5,
199 -1.13933444367982507207E 199 -1.13933444367982507207E6, -2.53252307177582951285E6, -2.01889141433532773231E6 };
200 static double LS2PI = 0.91893853320467274178 200 static double LS2PI = 0.91893853320467274178; /* log( sqrt( 2*pi ) ) */
201 #define MAXLGM 2.556348e305 201 #define MAXLGM 2.556348e305
202 202
203 /* 203 /*
204 ********************************************** 204 ************************************************************
205 */ 205 */
206 double nf_logGammaFunction( double x, nfu_stat 206 double nf_logGammaFunction( double x, nfu_status *status ) {
207 /* Logarithm of gamma function */ 207 /* Logarithm of gamma function */
208 208
209 int sgngam; 209 int sgngam;
210 210
211 *status = nfu_badInput; 211 *status = nfu_badInput;
212 if( !isfinite( x ) ) return( x ); 212 if( !isfinite( x ) ) return( x );
213 *status = nfu_Okay; 213 *status = nfu_Okay;
214 return( lgam( x, &sgngam, status ) ); 214 return( lgam( x, &sgngam, status ) );
215 } 215 }
216 /* 216 /*
217 ********************************************** 217 ************************************************************
218 */ 218 */
219 static double lgam( double x, int *sgngam, nfu 219 static double lgam( double x, int *sgngam, nfu_status *status ) {
220 220
221 double p, q, u, w, z; 221 double p, q, u, w, z;
222 int i; 222 int i;
223 223
224 *sgngam = 1; 224 *sgngam = 1;
225 225
226 if( x < -34.0 ) { 226 if( x < -34.0 ) {
227 q = -x; 227 q = -x;
228 w = lgam( q, sgngam, status ); 228 w = lgam( q, sgngam, status ); /* note this modifies *sgngam! */
229 p = floor( q ); 229 p = floor( q );
230 if( p == q ) goto lgsing; 230 if( p == q ) goto lgsing;
231 i = (int) p; 231 i = (int) p;
232 if( ( i & 1 ) == 0 ) { 232 if( ( i & 1 ) == 0 ) {
233 *sgngam = -1; } 233 *sgngam = -1; }
234 else { 234 else {
235 *sgngam = 1; 235 *sgngam = 1;
236 } 236 }
237 z = q - p; 237 z = q - p;
238 if( z > 0.5 ) { 238 if( z > 0.5 ) {
239 p += 1.0; 239 p += 1.0;
240 z = p - q; 240 z = p - q;
241 } 241 }
242 z = q * sin( M_PI * z ); 242 z = q * sin( M_PI * z );
243 if( z == 0.0 ) goto lgsing; 243 if( z == 0.0 ) goto lgsing;
244 z = LOGPI - G4Log( z ) - w; 244 z = LOGPI - G4Log( z ) - w;
245 return( z ); 245 return( z );
246 } 246 }
247 247
248 if( x < 13.0 ) { 248 if( x < 13.0 ) {
249 z = 1.0; 249 z = 1.0;
250 p = 0.0; 250 p = 0.0;
251 u = x; 251 u = x;
252 while( u >= 3.0 ) { 252 while( u >= 3.0 ) {
253 p -= 1.0; 253 p -= 1.0;
254 u = x + p; 254 u = x + p;
255 z *= u; 255 z *= u;
256 } // Loop checking, 11.06.2015, T. Koi 256 } // Loop checking, 11.06.2015, T. Koi
257 while( u < 2.0 ) { 257 while( u < 2.0 ) {
258 if( u == 0.0 ) goto lgsing; 258 if( u == 0.0 ) goto lgsing;
259 z /= u; 259 z /= u;
260 p += 1.0; 260 p += 1.0;
261 u = x + p; 261 u = x + p;
262 } // Loop checking, 11.06.2015, T. Koi 262 } // Loop checking, 11.06.2015, T. Koi
263 if( z < 0.0 ) { 263 if( z < 0.0 ) {
264 *sgngam = -1; 264 *sgngam = -1;
265 z = -z; } 265 z = -z; }
266 else { 266 else {
267 *sgngam = 1; 267 *sgngam = 1;
268 } 268 }
269 if( u == 2.0 ) return( G4Log( z ) ); 269 if( u == 2.0 ) return( G4Log( z ) );
270 p -= 2.0; 270 p -= 2.0;
271 x = x + p; 271 x = x + p;
272 p = x * nf_polevl( x, B, 5 ) / nf_p1ev 272 p = x * nf_polevl( x, B, 5 ) / nf_p1evl( x, C, 6);
273 return( G4Log( z ) + p ); 273 return( G4Log( z ) + p );
274 } 274 }
275 275
276 if( x > MAXLGM ) goto lgsing; 276 if( x > MAXLGM ) goto lgsing;
277 q = ( x - 0.5 ) * G4Log( x ) - x + LS2PI; 277 q = ( x - 0.5 ) * G4Log( x ) - x + LS2PI;
278 if( x > 1.0e8 ) return( q ); 278 if( x > 1.0e8 ) return( q );
279 279
280 p = 1.0 / ( x * x ); 280 p = 1.0 / ( x * x );
281 if( x >= 1000.0 ) { 281 if( x >= 1000.0 ) {
282 q += ( ( 7.9365079365079365079365e-4 * 282 q += ( ( 7.9365079365079365079365e-4 * p - 2.7777777777777777777778e-3 ) * p + 0.0833333333333333333333 ) / x; }
283 else { 283 else {
284 q += nf_polevl( p, A, 4 ) / x; 284 q += nf_polevl( p, A, 4 ) / x;
285 } 285 }
286 return( q ); 286 return( q );
287 287
288 lgsing: 288 lgsing:
289 return( *sgngam * DBL_MAX ); 289 return( *sgngam * DBL_MAX );
290 } 290 }
291 291
292 #if defined __cplusplus 292 #if defined __cplusplus
293 } 293 }
294 #endif 294 #endif
295 295