diff options
Diffstat (limited to 'libm/float/igamf.c')
| -rw-r--r-- | libm/float/igamf.c | 223 |
1 files changed, 0 insertions, 223 deletions
diff --git a/libm/float/igamf.c b/libm/float/igamf.c deleted file mode 100644 index c54225df4..000000000 --- a/libm/float/igamf.c +++ /dev/null @@ -1,223 +0,0 @@ -/* igamf.c - * - * Incomplete gamma integral - * - * - * - * SYNOPSIS: - * - * float a, x, y, igamf(); - * - * y = igamf( a, x ); - * - * - * - * DESCRIPTION: - * - * The function is defined by - * - * x - * - - * 1 | | -t a-1 - * igam(a,x) = ----- | e t dt. - * - | | - * | (a) - - * 0 - * - * - * In this implementation both arguments must be positive. - * The integral is evaluated by either a power series or - * continued fraction expansion, depending on the relative - * values of a and x. - * - * - * - * ACCURACY: - * - * Relative error: - * arithmetic domain # trials peak rms - * IEEE 0,30 20000 7.8e-6 5.9e-7 - * - */ - /* igamcf() - * - * Complemented incomplete gamma integral - * - * - * - * SYNOPSIS: - * - * float a, x, y, igamcf(); - * - * y = igamcf( a, x ); - * - * - * - * DESCRIPTION: - * - * The function is defined by - * - * - * igamc(a,x) = 1 - igam(a,x) - * - * inf. - * - - * 1 | | -t a-1 - * = ----- | e t dt. - * - | | - * | (a) - - * x - * - * - * In this implementation both arguments must be positive. - * The integral is evaluated by either a power series or - * continued fraction expansion, depending on the relative - * values of a and x. - * - * - * - * ACCURACY: - * - * Relative error: - * arithmetic domain # trials peak rms - * IEEE 0,30 30000 7.8e-6 5.9e-7 - * - */ - -/* -Cephes Math Library Release 2.2: June, 1992 -Copyright 1985, 1987, 1992 by Stephen L. Moshier -Direct inquiries to 30 Frost Street, Cambridge, MA 02140 -*/ - -#include <math.h> - -/* BIG = 1/MACHEPF */ -#define BIG 16777216. - -extern float MACHEPF, MAXLOGF; - -#ifdef ANSIC -float lgamf(float), expf(float), logf(float), igamf(float, float); -#else -float lgamf(), expf(), logf(), igamf(); -#endif - -#define fabsf(x) ( (x) < 0 ? -(x) : (x) ) - - - -float igamcf( float aa, float xx ) -{ -float a, x, ans, c, yc, ax, y, z; -float pk, pkm1, pkm2, qk, qkm1, qkm2; -float r, t; -static float big = BIG; - -a = aa; -x = xx; -if( (x <= 0) || ( a <= 0) ) - return( 1.0 ); - -if( (x < 1.0) || (x < a) ) - return( 1.0 - igamf(a,x) ); - -ax = a * logf(x) - x - lgamf(a); -if( ax < -MAXLOGF ) - { - mtherr( "igamcf", UNDERFLOW ); - return( 0.0 ); - } -ax = expf(ax); - -/* continued fraction */ -y = 1.0 - a; -z = x + y + 1.0; -c = 0.0; -pkm2 = 1.0; -qkm2 = x; -pkm1 = x + 1.0; -qkm1 = z * x; -ans = pkm1/qkm1; - -do - { - c += 1.0; - y += 1.0; - z += 2.0; - yc = y * c; - pk = pkm1 * z - pkm2 * yc; - qk = qkm1 * z - qkm2 * yc; - if( qk != 0 ) - { - r = pk/qk; - t = fabsf( (ans - r)/r ); - ans = r; - } - else - t = 1.0; - pkm2 = pkm1; - pkm1 = pk; - qkm2 = qkm1; - qkm1 = qk; - if( fabsf(pk) > big ) - { - pkm2 *= MACHEPF; - pkm1 *= MACHEPF; - qkm2 *= MACHEPF; - qkm1 *= MACHEPF; - } - } -while( t > MACHEPF ); - -return( ans * ax ); -} - - - -/* left tail of incomplete gamma function: - * - * inf. k - * a -x - x - * x e > ---------- - * - - - * k=0 | (a+k+1) - * - */ - -float igamf( float aa, float xx ) -{ -float a, x, ans, ax, c, r; - -a = aa; -x = xx; -if( (x <= 0) || ( a <= 0) ) - return( 0.0 ); - -if( (x > 1.0) && (x > a ) ) - return( 1.0 - igamcf(a,x) ); - -/* Compute x**a * exp(-x) / gamma(a) */ -ax = a * logf(x) - x - lgamf(a); -if( ax < -MAXLOGF ) - { - mtherr( "igamf", UNDERFLOW ); - return( 0.0 ); - } -ax = expf(ax); - -/* power series */ -r = a; -c = 1.0; -ans = 1.0; - -do - { - r += 1.0; - c *= x/r; - ans += c; - } -while( c/ans > MACHEPF ); - -return( ans * ax/a ); -} |