summaryrefslogtreecommitdiff
path: root/libm/ldouble/igaml.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
commit7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch)
tree3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/ldouble/igaml.c
parentc117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff)
Totally rework the math library, this time based on the MacOs X
math library (which is itself based on the math lib from FreeBSD). -Erik
Diffstat (limited to 'libm/ldouble/igaml.c')
-rw-r--r--libm/ldouble/igaml.c220
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 );
-}