diff options
Diffstat (limited to 'libm/ldouble/igaml.c')
| -rw-r--r-- | libm/ldouble/igaml.c | 220 |
1 files changed, 220 insertions, 0 deletions
diff --git a/libm/ldouble/igaml.c b/libm/ldouble/igaml.c new file mode 100644 index 000000000..0e59c5404 --- /dev/null +++ b/libm/ldouble/igaml.c @@ -0,0 +1,220 @@ +/* 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 ); +} |