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

1 files changed, 0 insertions, 223 deletions

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