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