Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
nf_gammaFunctions.cc File Reference
Include dependency graph for nf_gammaFunctions.cc:

Go to the source code of this file.

Macros

#define MAXGAM   171.624376956302725
 
#define MAXSTIR   143.01608
 
#define MAXLGM   2.556348e305
 

Functions

static double stirf (double x, nfu_status *status)
 
static double lgam (double x, int *sgngam, nfu_status *status)
 
double nf_gammaFunction (double x, nfu_status *status)
 
double nf_logGammaFunction (double x, nfu_status *status)
 

Variables

static double P []
 
static double Q []
 
static double LOGPI = 1.14472988584940017414
 
static double SQTPI = 2.50662827463100050242E0
 
static double STIR [5] = { 7.873113957930936284e-4, -2.2954996161337812638e-4, -2.6813261780578123283e-3, 3.472222216054586673e-3, 8.3333333333348225713e-2 }
 
static double A []
 
static double B []
 
static double C []
 
static double LS2PI = 0.91893853320467274178
 

Macro Definition Documentation

#define MAXGAM   171.624376956302725

Definition at line 93 of file nf_gammaFunctions.cc.

#define MAXLGM   2.556348e305

Definition at line 201 of file nf_gammaFunctions.cc.

#define MAXSTIR   143.01608

Definition at line 99 of file nf_gammaFunctions.cc.

Function Documentation

static double lgam ( double  x,
int sgngam,
nfu_status status 
)
static

Definition at line 219 of file nf_gammaFunctions.cc.

219  {
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 }
#define MAXLGM
#define M_PI
Definition: SbMath.h:34
const char * p
Definition: xmltok.h:285
double nf_polevl(double x, double coef[], int N)
Definition: nf_polevl.cc:46
tuple x
Definition: test.py:50
static double C[]
static double LOGPI
static double B[]
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
static double LS2PI
double nf_p1evl(double x, double coef[], int N)
Definition: nf_polevl.cc:67
static const G4double A[nN]
G4double G4Log(G4double x)
Definition: G4Log.hh:230
static double lgam(double x, int *sgngam, nfu_status *status)
tuple z
Definition: test.py:28
#define DBL_MAX
Definition: templates.hh:83

Here is the call graph for this function:

Here is the caller graph for this function:

double nf_gammaFunction ( double  x,
nfu_status status 
)

Definition at line 126 of file nf_gammaFunctions.cc.

126  {
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 }
#define M_PI
Definition: SbMath.h:34
const char * p
Definition: xmltok.h:285
double nf_polevl(double x, double coef[], int N)
Definition: nf_polevl.cc:46
static double Q[]
tuple x
Definition: test.py:50
static double stirf(double x, nfu_status *status)
static double P[]
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
#define isfinite
tuple z
Definition: test.py:28
#define DBL_MAX
Definition: templates.hh:83

Here is the call graph for this function:

Here is the caller graph for this function:

double nf_logGammaFunction ( double  x,
nfu_status status 
)

Definition at line 206 of file nf_gammaFunctions.cc.

206  {
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 }
tuple x
Definition: test.py:50
#define isfinite
static double lgam(double x, int *sgngam, nfu_status *status)

Here is the call graph for this function:

static double stirf ( double  x,
nfu_status status 
)
static

Definition at line 106 of file nf_gammaFunctions.cc.

106  {
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 }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
#define MAXSTIR
double nf_polevl(double x, double coef[], int N)
Definition: nf_polevl.cc:46
tuple x
Definition: test.py:50
static double SQTPI
static double STIR[5]
tuple v
Definition: test.py:18
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

double A[]
static
Initial value:
= { 8.11614167470508450300E-4, -5.95061904284301438324E-4, 7.93650340457716943945E-4,
-2.77777777730099687205E-3, 8.33333333333331927722E-2 }

Definition at line 194 of file nf_gammaFunctions.cc.

double B[]
static
Initial value:
= { -1.37825152569120859100E3, -3.88016315134637840924E4, -3.31612992738871184744E5,
-1.16237097492762307383E6, -1.72173700820839662146E6, -8.53555664245765465627E5 }

Definition at line 196 of file nf_gammaFunctions.cc.

double C[]
static
Initial value:
= { -3.51815701436523470549E2, -1.70642106651881159223E4, -2.20528590553854454839E5,
-1.13933444367982507207E6, -2.53252307177582951285E6, -2.01889141433532773231E6 }

Definition at line 198 of file nf_gammaFunctions.cc.

double LOGPI = 1.14472988584940017414
static

Definition at line 94 of file nf_gammaFunctions.cc.

double LS2PI = 0.91893853320467274178
static

Definition at line 200 of file nf_gammaFunctions.cc.

double P[]
static
Initial value:
= { 1.60119522476751861407E-4, 1.19135147006586384913E-3, 1.04213797561761569935E-2, 4.76367800457137231464E-2,
2.07448227648435975150E-1, 4.94214826801497100753E-1, 9.99999999999999996796E-1 }

Definition at line 89 of file nf_gammaFunctions.cc.

double Q[]
static
Initial value:
= { -2.31581873324120129819E-5, 5.39605580493303397842E-4, -4.45641913851797240494E-3, 1.18139785222060435552E-2,
3.58236398605498653373E-2, -2.34591795718243348568E-1, 7.14304917030273074085E-2, 1.00000000000000000320E0 }

Definition at line 91 of file nf_gammaFunctions.cc.

double SQTPI = 2.50662827463100050242E0
static

Definition at line 95 of file nf_gammaFunctions.cc.

double STIR[5] = { 7.873113957930936284e-4, -2.2954996161337812638e-4, -2.6813261780578123283e-3, 3.472222216054586673e-3, 8.3333333333348225713e-2 }
static

Definition at line 98 of file nf_gammaFunctions.cc.