summaryrefslogtreecommitdiff
path: root/libm/ldouble/igaml.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
committerEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
commit1077fa4d772832f77a677ce7fb7c2d513b959e3f (patch)
tree579bee13fb0b58d2800206366ec2caecbb15f3fc /libm/ldouble/igaml.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/ldouble/igaml.c')
-rw-r--r--libm/ldouble/igaml.c220
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 );
+}