78 #include "nf_specialFunctions.h"
80 #if defined __cplusplus
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;
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
101 static double stirf(
double x, nfu_status *status );
102 static double lgam(
double x,
int *sgngam, nfu_status *status );
106 static double stirf(
double x, nfu_status * ) {
131 *status = nfu_badInput;
140 if( p == q )
goto goverf;
142 if( ( i & 1 ) == 0 ) sgngam = -1;
148 z = q * sin( M_PI * z );
149 if( z == 0.0 )
goto goverf;
150 z = M_PI / ( fabs( z ) *
stirf( q, status ) );
153 z =
stirf( x, status );
155 return( sgngam * z );
165 if( x > -1.E-9 )
goto small;
171 if( x < 1.e-9 )
goto small;
176 if( x == 2.0 )
return( z );
184 if( x == 0.0 )
goto goverf;
185 return( z / ( ( 1.0 + 0.5772156649015329 * x ) * x ) );
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;
201 #define MAXLGM 2.556348e305
211 *status = nfu_badInput;
214 return(
lgam( x, &sgngam, status ) );
219 static double lgam(
double x,
int *sgngam, nfu_status *status ) {
221 double p, q, u,
w,
z;
228 w =
lgam( q, sgngam, status );
230 if( p == q )
goto lgsing;
232 if( ( i & 1 ) == 0 ) {
242 z = q * sin( M_PI * z );
243 if( z == 0.0 )
goto lgsing;
244 z = LOGPI -
G4Log( z ) -
w;
258 if( u == 0.0 )
goto lgsing;
269 if( u == 2.0 )
return(
G4Log( z ) );
273 return(
G4Log( z ) + p );
276 if( x >
MAXLGM )
goto lgsing;
278 if( x > 1.0e8 )
return( q );
282 q += ( ( 7.9365079365079365079365e-4 * p - 2.7777777777777777777778e-3 ) * p + 0.0833333333333333333333 ) /
x; }
292 #if defined __cplusplus
static G4Pow * GetInstance()
G4double powA(G4double A, G4double y) const
const G4double w[NPOINTSGL]
static double stirf(double x, nfu_status *status)
double nf_polevl(double x, double coef[], int N)
double nf_gammaFunction(double x, nfu_status *status)
G4double G4Log(G4double x)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
static double lgam(double x, int *sgngam, nfu_status *status)
const G4double x[NPOINTSGL]
double nf_logGammaFunction(double x, nfu_status *status)
double nf_p1evl(double x, double coef[], int N)