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

diff --git a/libm/float/hyp2f1f.c b/libm/float/hyp2f1f.c
deleted file mode 100644
index 01fe54928..000000000
--- a/libm/float/hyp2f1f.c
+++ /dev/null
@@ -1,442 +0,0 @@
-/*							hyp2f1f.c
- *
- *	Gauss hypergeometric function   F
- *	                               2 1
- *
- *
- * SYNOPSIS:
- *
- * float a, b, c, x, y, hyp2f1f();
- *
- * y = hyp2f1f( a, b, c, x );
- *
- *
- * DESCRIPTION:
- *
- *
- *  hyp2f1( a, b, c, x )  =   F ( a, b; c; x )
- *                           2 1
- *
- *           inf.
- *            -   a(a+1)...(a+k) b(b+1)...(b+k)   k+1
- *   =  1 +   >   -----------------------------  x   .
- *            -         c(c+1)...(c+k) (k+1)!
- *          k = 0
- *
- *  Cases addressed are
- *	Tests and escapes for negative integer a, b, or c
- *	Linear transformation if c - a or c - b negative integer
- *	Special case c = a or c = b
- *	Linear transformation for  x near +1
- *	Transformation for x < -0.5
- *	Psi function expansion if x > 0.5 and c - a - b integer
- *      Conditionally, a recurrence on c to make c-a-b > 0
- *
- * |x| > 1 is rejected.
- *
- * The parameters a, b, c are considered to be integer
- * valued if they are within 1.0e-6 of the nearest integer.
- *
- * ACCURACY:
- *
- *                      Relative error (-1 < x < 1):
- * arithmetic   domain     # trials      peak         rms
- *    IEEE      0,3         30000       5.8e-4      4.3e-6
- */
-
-/*							hyp2f1	*/
-
-
-/*
-Cephes Math Library Release 2.2:  November, 1992
-Copyright 1984, 1987, 1992 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
-*/
-
-
-#include <math.h>
-
-#define EPS 1.0e-5
-#define EPS2 1.0e-5
-#define ETHRESH 1.0e-5
-
-extern float MAXNUMF, MACHEPF;
-
-#define fabsf(x) ( (x) < 0 ? -(x) : (x) )
-
-#ifdef ANSIC
-float powf(float, float);
-static float hys2f1f(float, float, float, float, float *);
-static float hyt2f1f(float, float, float, float, float *);
-float gammaf(float), logf(float), expf(float), psif(float);
-float floorf(float);
-#else
-float powf(), gammaf(), logf(), expf(), psif();
-float floorf();
-static float hyt2f1f(), hys2f1f();
-#endif
-
-#define roundf(x) (floorf((x)+(float )0.5))
-
-
-
-
-float hyp2f1f( float aa, float bb, float cc, float xx )
-{
-float a, b, c, x;
-float d, d1, d2, e;
-float p, q, r, s, y, ax;
-float ia, ib, ic, id, err;
-int flag, i, aid;
-
-a = aa;
-b = bb;
-c = cc;
-x = xx;
-err = 0.0;
-ax = fabsf(x);
-s = 1.0 - x;
-flag = 0;
-ia = roundf(a); /* nearest integer to a */
-ib = roundf(b);
-
-if( a <= 0 )
-	{
-	if( fabsf(a-ia) < EPS )		/* a is a negative integer */
-		flag |= 1;
-	}
-
-if( b <= 0 )
-	{
-	if( fabsf(b-ib) < EPS )		/* b is a negative integer */
-		flag |= 2;
-	}
-
-if( ax < 1.0 )
-	{
-	if( fabsf(b-c) < EPS )		/* b = c */
-		{
-		y = powf( s, -a );	/* s to the -a power */
-		goto hypdon;
-		}
-	if( fabsf(a-c) < EPS )		/* a = c */
-		{
-		y = powf( s, -b );	/* s to the -b power */
-		goto hypdon;
-		}
-	}
-
-
-
-if( c <= 0.0 )
-	{
-	ic = roundf(c); 	/* nearest integer to c */
-	if( fabsf(c-ic) < EPS )		/* c is a negative integer */
-		{
-		/* check if termination before explosion */
-		if( (flag & 1) && (ia > ic) )
-			goto hypok;
-		if( (flag & 2) && (ib > ic) )
-			goto hypok;
-		goto hypdiv;
-		}
-	}
-
-if( flag )			/* function is a polynomial */
-	goto hypok;
-
-if( ax > 1.0 )			/* series diverges	*/
-	goto hypdiv;
-
-p = c - a;
-ia = roundf(p);
-if( (ia <= 0.0) && (fabsf(p-ia) < EPS) )	/* negative int c - a */
-	flag |= 4;
-
-r = c - b;
-ib = roundf(r); /* nearest integer to r */
-if( (ib <= 0.0) && (fabsf(r-ib) < EPS) )	/* negative int c - b */
-	flag |= 8;
-
-d = c - a - b;
-id = roundf(d); /* nearest integer to d */
-q = fabsf(d-id);
-
-if( fabsf(ax-1.0) < EPS )			/* |x| == 1.0	*/
-	{
-	if( x > 0.0 )
-		{
-		if( flag & 12 ) /* negative int c-a or c-b */
-			{
-			if( d >= 0.0 )
-				goto hypf;
-			else
-				goto hypdiv;
-			}
-		if( d <= 0.0 )
-			goto hypdiv;
-		y = gammaf(c)*gammaf(d)/(gammaf(p)*gammaf(r));
-		goto hypdon;
-		}
-
-	if( d <= -1.0 )
-		goto hypdiv;
-	}
-
-/* Conditionally make d > 0 by recurrence on c
- * AMS55 #15.2.27
- */
-if( d < 0.0 )
-	{
-/* Try the power series first */
-	y = hyt2f1f( a, b, c, x, &err );
-	if( err < ETHRESH )
-		goto hypdon;
-/* Apply the recurrence if power series fails */
-	err = 0.0;
-	aid = 2 - id;
-	e = c + aid;
-	d2 = hyp2f1f(a,b,e,x);
-	d1 = hyp2f1f(a,b,e+1.0,x);
-	q = a + b + 1.0;
-	for( i=0; i<aid; i++ )
-		{
-		r = e - 1.0;
-		y = (e*(r-(2.0*e-q)*x)*d2 + (e-a)*(e-b)*x*d1)/(e*r*s);
-		e = r;
-		d1 = d2;
-		d2 = y;
-		}
-	goto hypdon;
-	}
-
-
-if( flag & 12 )
-	goto hypf; /* negative integer c-a or c-b */
-
-hypok:
-y = hyt2f1f( a, b, c, x, &err );
-
-hypdon:
-if( err > ETHRESH )
-	{
-	mtherr( "hyp2f1", PLOSS );
-/*	printf( "Estimated err = %.2e\n", err );*/
-	}
-return(y);
-
-/* The transformation for c-a or c-b negative integer
- * AMS55 #15.3.3
- */
-hypf:
-y = powf( s, d ) * hys2f1f( c-a, c-b, c, x, &err );
-goto hypdon;
-
-/* The alarm exit */
-hypdiv:
-mtherr( "hyp2f1f", OVERFLOW );
-return( MAXNUMF );
-}
-
-
-
-
-/* Apply transformations for |x| near 1
- * then call the power series
- */
-static float hyt2f1f( float aa, float bb, float cc, float xx, float *loss )
-{
-float a, b, c, x;
-float p, q, r, s, t, y, d, err, err1;
-float ax, id, d1, d2, e, y1;
-int i, aid;
-
-a = aa;
-b = bb;
-c = cc;
-x = xx;
-err = 0.0;
-s = 1.0 - x;
-if( x < -0.5 )
-	{
-	if( b > a )
-		y = powf( s, -a ) * hys2f1f( a, c-b, c, -x/s, &err );
-
-	else
-		y = powf( s, -b ) * hys2f1f( c-a, b, c, -x/s, &err );
-
-	goto done;
-	}
-
-
-
-d = c - a - b;
-id = roundf(d);	/* nearest integer to d */
-
-if( x > 0.8 )
-{
-
-if( fabsf(d-id) > EPS2 ) /* test for integer c-a-b */
-	{
-/* Try the power series first */
-	y = hys2f1f( a, b, c, x, &err );
-	if( err < ETHRESH )
-		goto done;
-/* If power series fails, then apply AMS55 #15.3.6 */
-	q = hys2f1f( a, b, 1.0-d, s, &err );	
-	q *= gammaf(d) /(gammaf(c-a) * gammaf(c-b));
-	r = powf(s,d) * hys2f1f( c-a, c-b, d+1.0, s, &err1 );
-	r *= gammaf(-d)/(gammaf(a) * gammaf(b));
-	y = q + r;
-
-	q = fabsf(q); /* estimate cancellation error */
-	r = fabsf(r);
-	if( q > r )
-		r = q;
-	err += err1 + (MACHEPF*r)/y;
-
-	y *= gammaf(c);
-	goto done;
-	}	
-else
-	{
-/* Psi function expansion, AMS55 #15.3.10, #15.3.11, #15.3.12 */
-	if( id >= 0.0 )
-		{
-		e = d;
-		d1 = d;
-		d2 = 0.0;
-		aid = id;
-		}
-	else
-		{
-		e = -d;
-		d1 = 0.0;
-		d2 = d;
-		aid = -id;
-		}
-
-	ax = logf(s);
-
-	/* sum for t = 0 */
-	y = psif(1.0) + psif(1.0+e) - psif(a+d1) - psif(b+d1) - ax;
-	y /= gammaf(e+1.0);
-
-	p = (a+d1) * (b+d1) * s / gammaf(e+2.0);	/* Poch for t=1 */
-	t = 1.0;
-	do
-		{
-		r = psif(1.0+t) + psif(1.0+t+e) - psif(a+t+d1)
-			- psif(b+t+d1) - ax;
-		q = p * r;
-		y += q;
-		p *= s * (a+t+d1) / (t+1.0);
-		p *= (b+t+d1) / (t+1.0+e);
-		t += 1.0;
-		}
-	while( fabsf(q/y) > EPS );
-
-
-	if( id == 0.0 )
-		{
-		y *= gammaf(c)/(gammaf(a)*gammaf(b));
-		goto psidon;
-		}
-
-	y1 = 1.0;
-
-	if( aid == 1 )
-		goto nosum;
-
-	t = 0.0;
-	p = 1.0;
-	for( i=1; i<aid; i++ )
-		{
-		r = 1.0-e+t;
-		p *= s * (a+t+d2) * (b+t+d2) / r;
-		t += 1.0;
-		p /= t;
-		y1 += p;
-		}
-
-
-nosum:
-	p = gammaf(c);
-	y1 *= gammaf(e) * p / (gammaf(a+d1) * gammaf(b+d1));
-	y *= p / (gammaf(a+d2) * gammaf(b+d2));
-	if( (aid & 1) != 0 )
-		y = -y;
-
-	q = powf( s, id );	/* s to the id power */
-	if( id > 0.0 )
-		y *= q;
-	else
-		y1 *= q;
-
-	y += y1;
-psidon:
-	goto done;
-	}
-}
-
-
-/* Use defining power series if no special cases */
-y = hys2f1f( a, b, c, x, &err );
-
-done:
-*loss = err;
-return(y);
-}
-
-
-
-
-
-/* Defining power series expansion of Gauss hypergeometric function */
-
-static float hys2f1f( float aa, float bb, float cc, float xx, float *loss )
-{
-int i;
-float a, b, c, x;
-float f, g, h, k, m, s, u, umax;
-
-
-a = aa;
-b = bb;
-c = cc;
-x = xx;
-i = 0;
-umax = 0.0;
-f = a;
-g = b;
-h = c;
-k = 0.0;
-s = 1.0;
-u = 1.0;
-
-do
-	{
-	if( fabsf(h) < EPS )
-		return( MAXNUMF );
-	m = k + 1.0;
-	u = u * ((f+k) * (g+k) * x / ((h+k) * m));
-	s += u;
-	k = fabsf(u);  /* remember largest term summed */
-	if( k > umax )
-		umax = k;
-	k = m;
-	if( ++i > 10000 ) /* should never happen */
-		{
-		*loss = 1.0;
-		return(s);
-		}
-	}
-while( fabsf(u/s) > MACHEPF );
-
-/* return estimated relative error */
-*loss = (MACHEPF*umax)/fabsf(s) + (MACHEPF*i);
-
-return(s);
-}
-
-