Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lend/src/nf_gammaFunctions.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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
295