uClibc now has a math library. muahahahaha!

 -Erik

1 files changed, 223 insertions, 0 deletions

diff --git a/libm/float/igamf.c b/libm/float/igamf.c
new file mode 100644
index 000000000..c54225df4
--- /dev/null
+++ b/libm/float/igamf.c
@@ -0,0 +1,223 @@
+/*							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 );
+}