diff options
Diffstat (limited to 'libm/ldouble/igaml.c')
| -rw-r--r-- | libm/ldouble/igaml.c | 220 |
1 files changed, 0 insertions, 220 deletions
diff --git a/libm/ldouble/igaml.c b/libm/ldouble/igaml.c deleted file mode 100644 index 0e59c5404..000000000 --- a/libm/ldouble/igaml.c +++ /dev/null @@ -1,220 +0,0 @@ -/* igaml.c - * - * Incomplete gamma integral - * - * - * - * SYNOPSIS: - * - * long double a, x, y, igaml(); - * - * y = igaml( 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 - * DEC 0,30 4000 4.4e-15 6.3e-16 - * IEEE 0,30 10000 3.6e-14 5.1e-15 - * - */ - /* igamcl() - * - * Complemented incomplete gamma integral - * - * - * - * SYNOPSIS: - * - * long double a, x, y, igamcl(); - * - * y = igamcl( 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 - * DEC 0,30 2000 2.7e-15 4.0e-16 - * IEEE 0,30 60000 1.4e-12 6.3e-15 - * - */ - -/* -Cephes Math Library Release 2.3: March, 1995 -Copyright 1985, 1995 by Stephen L. Moshier -*/ - -#include <math.h> -#ifdef ANSIPROT -extern long double lgaml ( long double ); -extern long double expl ( long double ); -extern long double logl ( long double ); -extern long double fabsl ( long double ); -extern long double gammal ( long double ); -long double igaml ( long double, long double ); -long double igamcl ( long double, long double ); -#else -long double lgaml(), expl(), logl(), fabsl(), igaml(), gammal(); -long double igamcl(); -#endif - -#define BIG 9.223372036854775808e18L -#define MAXGAML 1755.455L -extern long double MACHEPL, MINLOGL; - -long double igamcl( a, x ) -long double a, x; -{ -long double ans, c, yc, ax, y, z, r, t; -long double pk, pkm1, pkm2, qk, qkm1, qkm2; - -if( (x <= 0.0L) || ( a <= 0.0L) ) - return( 1.0L ); - -if( (x < 1.0L) || (x < a) ) - return( 1.0L - igaml(a,x) ); - -ax = a * logl(x) - x - lgaml(a); -if( ax < MINLOGL ) - { - mtherr( "igamcl", UNDERFLOW ); - return( 0.0L ); - } -ax = expl(ax); - -/* continued fraction */ -y = 1.0L - a; -z = x + y + 1.0L; -c = 0.0L; -pkm2 = 1.0L; -qkm2 = x; -pkm1 = x + 1.0L; -qkm1 = z * x; -ans = pkm1/qkm1; - -do - { - c += 1.0L; - y += 1.0L; - z += 2.0L; - yc = y * c; - pk = pkm1 * z - pkm2 * yc; - qk = qkm1 * z - qkm2 * yc; - if( qk != 0.0L ) - { - r = pk/qk; - t = fabsl( (ans - r)/r ); - ans = r; - } - else - t = 1.0L; - pkm2 = pkm1; - pkm1 = pk; - qkm2 = qkm1; - qkm1 = qk; - if( fabsl(pk) > BIG ) - { - pkm2 /= BIG; - pkm1 /= BIG; - qkm2 /= BIG; - qkm1 /= BIG; - } - } -while( t > MACHEPL ); - -return( ans * ax ); -} - - - -/* left tail of incomplete gamma function: - * - * inf. k - * a -x - x - * x e > ---------- - * - - - * k=0 | (a+k+1) - * - */ - -long double igaml( a, x ) -long double a, x; -{ -long double ans, ax, c, r; - -if( (x <= 0.0L) || ( a <= 0.0L) ) - return( 0.0L ); - -if( (x > 1.0L) && (x > a ) ) - return( 1.0L - igamcl(a,x) ); - -ax = a * logl(x) - x - lgaml(a); -if( ax < MINLOGL ) - { - mtherr( "igaml", UNDERFLOW ); - return( 0.0L ); - } -ax = expl(ax); - -/* power series */ -r = a; -c = 1.0L; -ans = 1.0L; - -do - { - r += 1.0L; - c *= x/r; - ans += c; - } -while( c/ans > MACHEPL ); - -return( ans * ax/a ); -} |