Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
nf_incompleteGammaFunctions.cc
Go to the documentation of this file.
1 /* igam.c
2  *
3  * Incomplete gamma integral
4  *
5  *
6  * DESCRIPTION:
7  *
8  * The function is defined by
9  *
10  * x
11  * -
12  * | | -t a-1
13  * igam(a,x) = | e t dt.
14  * | |
15  * -
16  * 0
17  *
18  *
19  * In this implementation both arguments must be positive.
20  * The integral is evaluated by either a power series or
21  * continued fraction expansion, depending on the relative
22  * values of a and x.
23  *
24  * ACCURACY:
25  *
26  * Relative error:
27  * arithmetic domain # trials peak rms
28  * IEEE 0,30 200000 3.6e-14 2.9e-15
29  * IEEE 0,100 300000 9.9e-14 1.5e-14
30  */
31 /* igamc()
32  *
33  * Complemented incomplete gamma integral
34  *
35  *
36  * DESCRIPTION:
37  *
38  * The function is defined by
39  *
40  *
41  * igamc(a,x) = 1 - igam(a,x)
42  *
43  * inf.
44  * -
45  * | | -t a-1
46  * = | e t dt.
47  * | |
48  * -
49  * x
50  *
51  *
52  * In this implementation both arguments must be positive.
53  * The integral is evaluated by either a power series or
54  * continued fraction expansion, depending on the relative
55  * values of a and x.
56  *
57  * ACCURACY:
58  *
59  * Tested at random a, x.
60  * a x Relative error:
61  * arithmetic domain domain # trials peak rms
62  * IEEE 0.5,100 0,100 200000 1.9e-14 1.7e-15
63  * IEEE 0.01,0.5 0,100 200000 1.4e-13 1.6e-15
64  */
65 
66 /*
67 Cephes Math Library Release 2.8: June, 2000
68 Copyright 1985, 1987, 2000 by Stephen L. Moshier
69 */
70 
71 #include "nf_specialFunctions.h"
72 
73 #if defined __cplusplus
74 #include <cmath>
75 #include "G4Exp.hh"
76 #include "G4Log.hh"
77 namespace GIDI {
78 using namespace GIDI;
79 using namespace std;
80 #endif
81 
82 static double big = 4.503599627370496e15;
83 static double biginv = 2.22044604925031308085e-16;
84 
85 /*
86 ************************************************************
87 */
88 double nf_incompleteGammaFunctionComplementary( double a, double x, nfu_status *status ) {
89 
90  double ans, ax, c, yc, r, t, y, z;
91  double pk, pkm1, pkm2, qk, qkm1, qkm2;
92 
93  *status = nfu_badInput;
94  if( !isfinite( x ) ) return( x );
95  *status = nfu_Okay;
96 
97  if( ( x <= 0 ) || ( a <= 0 ) ) return( 1.0 );
98  if( ( x < 1.0 ) || ( x < a ) ) return( nf_gammaFunction( a, status ) - nf_incompleteGammaFunction( a, x, status ) );
99 
100  ax = G4Exp( a * G4Log( x ) - x );
101  if( ax == 0. ) return( 0.0 );
102 
103  if( x < 10000. ) {
104  y = 1.0 - a; /* continued fraction */
105  z = x + y + 1.0;
106  c = 0.0;
107  pkm2 = 1.0;
108  qkm2 = x;
109  pkm1 = x + 1.0;
110  qkm1 = z * x;
111  ans = pkm1 / qkm1;
112 
113  do {
114  c += 1.0;
115  y += 1.0;
116  z += 2.0;
117  yc = y * c;
118  pk = pkm1 * z - pkm2 * yc;
119  qk = qkm1 * z - qkm2 * yc;
120  if( qk != 0 ) {
121  r = pk / qk;
122  t = fabs( ( ans - r ) / r );
123  ans = r; }
124  else {
125  t = 1.0;
126  }
127  pkm2 = pkm1;
128  pkm1 = pk;
129  qkm2 = qkm1;
130  qkm1 = qk;
131  if( fabs( pk ) > big ) {
132  pkm2 *= biginv;
133  pkm1 *= biginv;
134  qkm2 *= biginv;
135  qkm1 *= biginv;
136  }
137  } while( t > DBL_EPSILON ); } // Loop checking, 11.06.2015, T. Koi
138  else { /* Asymptotic expansion. */
139  y = 1. / x;
140  r = a;
141  c = 1.;
142  ans = 1.;
143  do {
144  a -= 1.;
145  c *= a * y;
146  ans += c;
147  } while( fabs( c ) > 100 * ans * DBL_EPSILON ); // Loop checking, 11.06.2015, T. Koi
148  }
149 
150  return( ans * ax );
151 }
152 /*
153 ************************************************************
154 */
155 double nf_incompleteGammaFunction( double a, double x, nfu_status *status ) {
156 /* left tail of incomplete gamma function:
157 *
158 * inf. k
159 * a -x - x
160 * x e > ----------
161 * - -
162 * k=0 | (a+k+1)
163 */
164  double ans, ax, c, r;
165 
166  *status = nfu_badInput;
167  if( !isfinite( x ) ) return( x );
168  *status = nfu_Okay;
169 
170  if( ( x <= 0 ) || ( a <= 0 ) ) return( 0.0 );
171  if( ( x > 1.0 ) && ( x > a ) ) return( nf_gammaFunction( a, status ) - nf_incompleteGammaFunctionComplementary( a, x, status ) );
172 
173  ax = G4Exp( a * G4Log( x ) - x ); /* Compute x**a * exp(-x) */
174  if( ax == 0. ) return( 0.0 );
175 
176  r = a; /* power series */
177  c = 1.0;
178  ans = 1.0;
179  do {
180  r += 1.0;
181  c *= x / r;
182  ans += c;
183  } while( c > ans * DBL_EPSILON ); // Loop checking, 11.06.2015, T. Koi
184 
185  return( ans * ax / a );
186 }
187 
188 #if defined __cplusplus
189 }
190 #endif
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
tuple x
Definition: test.py:50
static double biginv
enum nfu_status_e nfu_status
#define DBL_EPSILON
Definition: templates.hh:87
#define isfinite
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
double nf_incompleteGammaFunctionComplementary(double a, double x, nfu_status *status)
tuple z
Definition: test.py:28
double nf_incompleteGammaFunction(double a, double x, nfu_status *status)
tuple c
Definition: test.py:13
static double big
double nf_gammaFunction(double x, nfu_status *status)