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, 424 deletions

diff --git a/libm/float/incbetf.c b/libm/float/incbetf.c
deleted file mode 100644
index fed9aae4b..000000000
--- a/libm/float/incbetf.c
+++ /dev/null
@@ -1,424 +0,0 @@
-/*							incbetf.c
- *
- *	Incomplete beta integral
- *
- *
- * SYNOPSIS:
- *
- * float a, b, x, y, incbetf();
- *
- * y = incbetf( a, b, x );
- *
- *
- * DESCRIPTION:
- *
- * Returns incomplete beta integral of the arguments, evaluated
- * from zero to x.  The function is defined as
- *
- *                  x
- *     -            -
- *    | (a+b)      | |  a-1     b-1
- *  -----------    |   t   (1-t)   dt.
- *   -     -     | |
- *  | (a) | (b)   -
- *                 0
- *
- * The domain of definition is 0 <= x <= 1.  In this
- * implementation a and b are restricted to positive values.
- * The integral from x to 1 may be obtained by the symmetry
- * relation
- *
- *    1 - incbet( a, b, x )  =  incbet( b, a, 1-x ).
- *
- * The integral is evaluated by a continued fraction expansion.
- * If a < 1, the function calls itself recursively after a
- * transformation to increase a to a+1.
- *
- * ACCURACY:
- *
- * Tested at random points (a,b,x) with a and b in the indicated
- * interval and x between 0 and 1.
- *
- * arithmetic   domain     # trials      peak         rms
- * Relative error:
- *    IEEE       0,30       10000       3.7e-5      5.1e-6
- *    IEEE       0,100      10000       1.7e-4      2.5e-5
- * The useful domain for relative error is limited by underflow
- * of the single precision exponential function.
- * Absolute error:
- *    IEEE       0,30      100000       2.2e-5      9.6e-7
- *    IEEE       0,100      10000       6.5e-5      3.7e-6
- *
- * Larger errors may occur for extreme ratios of a and b.
- *
- * ERROR MESSAGES:
- *   message         condition      value returned
- * incbetf domain     x<0, x>1          0.0
- */
-
-
-/*
-Cephes Math Library, Release 2.2:  July, 1992
-Copyright 1984, 1987, 1992 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
-*/
-
-#include <math.h>
-
-#ifdef ANSIC
-float lgamf(float), expf(float), logf(float);
-static float incbdf(float, float, float);
-static float incbcff(float, float, float);
-float incbpsf(float, float, float);
-#else
-float lgamf(), expf(), logf();
-float incbpsf();
-static float incbcff(), incbdf();
-#endif
-
-#define fabsf(x) ( (x) < 0 ? -(x) : (x) )
-
-/* BIG = 1/MACHEPF */
-#define BIG   16777216.
-extern float MACHEPF, MAXLOGF;
-#define MINLOGF (-MAXLOGF)
-
-float incbetf( float aaa, float bbb, float xxx )
-{
-float aa, bb, xx, ans, a, b, t, x, onemx;
-int flag;
-
-aa = aaa;
-bb = bbb;
-xx = xxx;
-if( (xx <= 0.0) || ( xx >= 1.0) )
-	{
-	if( xx == 0.0 )
-		return(0.0);
-	if( xx == 1.0 )
-		return( 1.0 );
-	mtherr( "incbetf", DOMAIN );
-	return( 0.0 );
-	}
-
-onemx = 1.0 - xx;
-
-
-/* transformation for small aa */
-
-if( aa <= 1.0 )
-	{
-	ans = incbetf( aa+1.0, bb, xx );
-	t = aa*logf(xx) + bb*logf( 1.0-xx )
-		+ lgamf(aa+bb) - lgamf(aa+1.0) - lgamf(bb);
-	if( t > MINLOGF )
-		ans += expf(t);
-	return( ans );
-	}
-
-
-/* see if x is greater than the mean */
-
-if( xx > (aa/(aa+bb)) )
-	{
-	flag = 1;
-	a = bb;
-	b = aa;
-	t = xx;
-	x = onemx;
-	}
-else
-	{
-	flag = 0;
-	a = aa;
-	b = bb;
-	t = onemx;
-	x = xx;
-	}
-
-/* transformation for small aa */
-/*
-if( a <= 1.0 )
-	{
-	ans = a*logf(x) + b*logf( onemx )
-		+ lgamf(a+b) - lgamf(a+1.0) - lgamf(b);
- 	t = incbetf( a+1.0, b, x );
-	if( ans > MINLOGF )
-		t += expf(ans);
-	goto bdone;
-	}
-*/
-/* Choose expansion for optimal convergence */
-
-
-if( b > 10.0 )
-	{
-if( fabsf(b*x/a) < 0.3 )
-	{
-	t = incbpsf( a, b, x );
-	goto bdone;
-	}
-	}
-
-ans = x * (a+b-2.0)/(a-1.0);
-if( ans < 1.0 )
-	{
-	ans = incbcff( a, b, x );
-	t = b * logf( t );
-	}
-else
-	{
-	ans = incbdf( a, b, x );
-	t = (b-1.0) * logf(t);
-	}
-
-t += a*logf(x) + lgamf(a+b) - lgamf(a) - lgamf(b);
-t += logf( ans/a );
-
-if( t < MINLOGF )
-	{
-	t = 0.0;
-	if( flag == 0 )
-		{
-		mtherr( "incbetf", UNDERFLOW );
-		}
-	}
-else
-	{
-	t = expf(t);
-	}
-bdone:
-
-if( flag )
-	t = 1.0 - t;
-
-return( t );
-}
-
-/* Continued fraction expansion #1
- * for incomplete beta integral
- */
-
-static float incbcff( float aa, float bb, float xx )
-{
-float a, b, x, xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
-float k1, k2, k3, k4, k5, k6, k7, k8;
-float r, t, ans;
-static float big = BIG;
-int n;
-
-a = aa;
-b = bb;
-x = xx;
-k1 = a;
-k2 = a + b;
-k3 = a;
-k4 = a + 1.0;
-k5 = 1.0;
-k6 = b - 1.0;
-k7 = k4;
-k8 = a + 2.0;
-
-pkm2 = 0.0;
-qkm2 = 1.0;
-pkm1 = 1.0;
-qkm1 = 1.0;
-ans = 1.0;
-r = 0.0;
-n = 0;
-do
-	{
-	
-	xk = -( x * k1 * k2 )/( k3 * k4 );
-	pk = pkm1 +  pkm2 * xk;
-	qk = qkm1 +  qkm2 * xk;
-	pkm2 = pkm1;
-	pkm1 = pk;
-	qkm2 = qkm1;
-	qkm1 = qk;
-
-	xk = ( x * k5 * k6 )/( k7 * k8 );
-	pk = pkm1 +  pkm2 * xk;
-	qk = qkm1 +  qkm2 * xk;
-	pkm2 = pkm1;
-	pkm1 = pk;
-	qkm2 = qkm1;
-	qkm1 = qk;
-
-	if( qk != 0 )
-		r = pk/qk;
-	if( r != 0 )
-		{
-		t = fabsf( (ans - r)/r );
-		ans = r;
-		}
-	else
-		t = 1.0;
-
-	if( t < MACHEPF )
-		goto cdone;
-
-	k1 += 1.0;
-	k2 += 1.0;
-	k3 += 2.0;
-	k4 += 2.0;
-	k5 += 1.0;
-	k6 -= 1.0;
-	k7 += 2.0;
-	k8 += 2.0;
-
-	if( (fabsf(qk) + fabsf(pk)) > big )
-		{
-		pkm2 *= MACHEPF;
-		pkm1 *= MACHEPF;
-		qkm2 *= MACHEPF;
-		qkm1 *= MACHEPF;
-		}
-	if( (fabsf(qk) < MACHEPF) || (fabsf(pk) < MACHEPF) )
-		{
-		pkm2 *= big;
-		pkm1 *= big;
-		qkm2 *= big;
-		qkm1 *= big;
-		}
-	}
-while( ++n < 100 );
-
-cdone:
-return(ans);
-}
-
-
-/* Continued fraction expansion #2
- * for incomplete beta integral
- */
-
-static float incbdf( float aa, float bb, float xx )
-{
-float a, b, x, xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
-float k1, k2, k3, k4, k5, k6, k7, k8;
-float r, t, ans, z;
-static float big = BIG;
-int n;
-
-a = aa;
-b = bb;
-x = xx;
-k1 = a;
-k2 = b - 1.0;
-k3 = a;
-k4 = a + 1.0;
-k5 = 1.0;
-k6 = a + b;
-k7 = a + 1.0;;
-k8 = a + 2.0;
-
-pkm2 = 0.0;
-qkm2 = 1.0;
-pkm1 = 1.0;
-qkm1 = 1.0;
-z = x / (1.0-x);
-ans = 1.0;
-r = 0.0;
-n = 0;
-do
-	{
-	
-	xk = -( z * k1 * k2 )/( k3 * k4 );
-	pk = pkm1 +  pkm2 * xk;
-	qk = qkm1 +  qkm2 * xk;
-	pkm2 = pkm1;
-	pkm1 = pk;
-	qkm2 = qkm1;
-	qkm1 = qk;
-
-	xk = ( z * k5 * k6 )/( k7 * k8 );
-	pk = pkm1 +  pkm2 * xk;
-	qk = qkm1 +  qkm2 * xk;
-	pkm2 = pkm1;
-	pkm1 = pk;
-	qkm2 = qkm1;
-	qkm1 = qk;
-
-	if( qk != 0 )
-		r = pk/qk;
-	if( r != 0 )
-		{
-		t = fabsf( (ans - r)/r );
-		ans = r;
-		}
-	else
-		t = 1.0;
-
-	if( t < MACHEPF )
-		goto cdone;
-
-	k1 += 1.0;
-	k2 -= 1.0;
-	k3 += 2.0;
-	k4 += 2.0;
-	k5 += 1.0;
-	k6 += 1.0;
-	k7 += 2.0;
-	k8 += 2.0;
-
-	if( (fabsf(qk) + fabsf(pk)) > big )
-		{
-		pkm2 *= MACHEPF;
-		pkm1 *= MACHEPF;
-		qkm2 *= MACHEPF;
-		qkm1 *= MACHEPF;
-		}
-	if( (fabsf(qk) < MACHEPF) || (fabsf(pk) < MACHEPF) )
-		{
-		pkm2 *= big;
-		pkm1 *= big;
-		qkm2 *= big;
-		qkm1 *= big;
-		}
-	}
-while( ++n < 100 );
-
-cdone:
-return(ans);
-}
-
-
-/* power series */
-float incbpsf( float aa, float bb, float xx )
-{
-float a, b, x, t, u, y, s;
-
-a = aa;
-b = bb;
-x = xx;
-
-y = a * logf(x) + (b-1.0)*logf(1.0-x) - logf(a);
-y -= lgamf(a) + lgamf(b);
-y += lgamf(a+b);
-
-
-t = x / (1.0 - x);
-s = 0.0;
-u = 1.0;
-do
-	{
-	b -= 1.0;
-	if( b == 0.0 )
-		break;
-	a += 1.0;
-	u *= t*b/a;
-	s += u;
-	}
-while( fabsf(u) > MACHEPF );
-
-if( y < MINLOGF )
-	{
-	mtherr( "incbetf", UNDERFLOW );
-	s = 0.0;
-	}
-else
-	s = expf(y) * (1.0 + s);
-/*printf( "incbpsf: %.4e\n", s );*/
-return(s);
-}