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