3 * Incomplete gamma integral
9 * float a, x, y, igamf();
17 * The function is defined by
22 * igam(a,x) = ----- | e t dt.
28 * In this implementation both arguments must be positive.
29 * The integral is evaluated by either a power series or
30 * continued fraction expansion, depending on the relative
38 * arithmetic domain # trials peak rms
39 * IEEE 0,30 20000 7.8e-6 5.9e-7
44 * Complemented incomplete gamma integral
50 * float a, x, y, igamcf();
58 * The function is defined by
61 * igamc(a,x) = 1 - igam(a,x)
72 * In this implementation both arguments must be positive.
73 * The integral is evaluated by either a power series or
74 * continued fraction expansion, depending on the relative
82 * arithmetic domain # trials peak rms
83 * IEEE 0,30 30000 7.8e-6 5.9e-7
88 Cephes Math Library Release 2.2: June, 1992
89 Copyright 1985, 1987, 1992 by Stephen L. Moshier
90 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
98 extern float MACHEPF, MAXLOGF;
101 float lgamf(float), expf(float), logf(float), igamf(float, float);
103 float lgamf(), expf(), logf(), igamf();
106 #define fabsf(x) ( (x) < 0 ? -(x) : (x) )
110 float igamcf( float aa, float xx )
112 float a, x, ans, c, yc, ax, y, z;
113 float pk, pkm1, pkm2, qk, qkm1, qkm2;
115 static float big = BIG;
119 if( (x <= 0) || ( a <= 0) )
122 if( (x < 1.0) || (x < a) )
123 return( 1.0 - igamf(a,x) );
125 ax = a * logf(x) - x - lgamf(a);
128 mtherr( "igamcf", UNDERFLOW );
133 /* continued fraction */
149 pk = pkm1 * z - pkm2 * yc;
150 qk = qkm1 * z - qkm2 * yc;
154 t = fabsf( (ans - r)/r );
163 if( fabsf(pk) > big )
171 while( t > MACHEPF );
178 /* left tail of incomplete gamma function:
188 float igamf( float aa, float xx )
190 float a, x, ans, ax, c, r;
194 if( (x <= 0) || ( a <= 0) )
197 if( (x > 1.0) && (x > a ) )
198 return( 1.0 - igamcf(a,x) );
200 /* Compute x**a * exp(-x) / gamma(a) */
201 ax = a * logf(x) - x - lgamf(a);
204 mtherr( "igamf", UNDERFLOW );
220 while( c/ans > MACHEPF );
222 return( ans * ax/a );