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

diff --git a/libm/float/hypergf.c b/libm/float/hypergf.c
deleted file mode 100644
index 60d0eb4c5..000000000
--- a/libm/float/hypergf.c
+++ /dev/null
@@ -1,384 +0,0 @@
-/*							hypergf.c
- *
- *	Confluent hypergeometric function
- *
- *
- *
- * SYNOPSIS:
- *
- * float a, b, x, y, hypergf();
- *
- * y = hypergf( a, b, x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Computes the confluent hypergeometric function
- *
- *                          1           2
- *                       a x    a(a+1) x
- *   F ( a,b;x )  =  1 + ---- + --------- + ...
- *  1 1                  b 1!   b(b+1) 2!
- *
- * Many higher transcendental functions are special cases of
- * this power series.
- *
- * As is evident from the formula, b must not be a negative
- * integer or zero unless a is an integer with 0 >= a > b.
- *
- * The routine attempts both a direct summation of the series
- * and an asymptotic expansion.  In each case error due to
- * roundoff, cancellation, and nonconvergence is estimated.
- * The result with smaller estimated error is returned.
- *
- *
- *
- * ACCURACY:
- *
- * Tested at random points (a, b, x), all three variables
- * ranging from 0 to 30.
- *                      Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    IEEE      0,5         10000       6.6e-7      1.3e-7
- *    IEEE      0,30        30000       1.1e-5      6.5e-7
- *
- * Larger errors can be observed when b is near a negative
- * integer or zero.  Certain combinations of arguments yield
- * serious cancellation error in the power series summation
- * and also are not in the region of near convergence of the
- * asymptotic series.  An error message is printed if the
- * self-estimated relative error is greater than 1.0e-3.
- *
- */
-
-/*							hyperg.c */
-
-
-/*
-Cephes Math Library Release 2.1:  November, 1988
-Copyright 1984, 1987, 1988 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
-*/
-
-#include <math.h>
-
-extern float MAXNUMF, MACHEPF;
-
-#define fabsf(x) ( (x) < 0 ? -(x) : (x) )
-
-#ifdef ANSIC
-float expf(float);
-float hyp2f0f(float, float, float, int, float *);
-static float hy1f1af(float, float, float, float *);
-static float hy1f1pf(float, float, float, float *);
-float logf(float), gammaf(float), lgamf(float);
-#else
-float expf(), hyp2f0f();
-float logf(), gammaf(), lgamf();
-static float hy1f1pf(), hy1f1af();
-#endif
-
-float hypergf( float aa, float bb, float xx )
-{
-float a, b, x, asum, psum, acanc, pcanc, temp;
-
-
-a = aa;
-b = bb;
-x = xx;
-/* See if a Kummer transformation will help */
-temp = b - a;
-if( fabsf(temp) < 0.001 * fabsf(a) )
-	return( expf(x) * hypergf( temp, b, -x )  );
-
-psum = hy1f1pf( a, b, x, &pcanc );
-if( pcanc < 1.0e-6 )
-	goto done;
-
-
-/* try asymptotic series */
-
-asum = hy1f1af( a, b, x, &acanc );
-
-
-/* Pick the result with less estimated error */
-
-if( acanc < pcanc )
-	{
-	pcanc = acanc;
-	psum = asum;
-	}
-
-done:
-if( pcanc > 1.0e-3 )
-	mtherr( "hyperg", PLOSS );
-
-return( psum );
-}
-
-
-
-
-/* Power series summation for confluent hypergeometric function		*/
-
-
-static float hy1f1pf( float aa, float bb, float xx, float *err )
-{
-float a, b, x, n, a0, sum, t, u, temp;
-float an, bn, maxt, pcanc;
-
-a = aa;
-b = bb;
-x = xx;
-/* set up for power series summation */
-an = a;
-bn = b;
-a0 = 1.0;
-sum = 1.0;
-n = 1.0;
-t = 1.0;
-maxt = 0.0;
-
-
-while( t > MACHEPF )
-	{
-	if( bn == 0 )			/* check bn first since if both	*/
-		{
-		mtherr( "hypergf", SING );
-		return( MAXNUMF );	/* an and bn are zero it is	*/
-		}
-	if( an == 0 )			/* a singularity		*/
-		return( sum );
-	if( n > 200 )
-		goto pdone;
-	u = x * ( an / (bn * n) );
-
-	/* check for blowup */
-	temp = fabsf(u);
-	if( (temp > 1.0 ) && (maxt > (MAXNUMF/temp)) )
-		{
-		pcanc = 1.0;	/* estimate 100% error */
-		goto blowup;
-		}
-
-	a0 *= u;
-	sum += a0;
-	t = fabsf(a0);
-	if( t > maxt )
-		maxt = t;
-/*
-	if( (maxt/fabsf(sum)) > 1.0e17 )
-		{
-		pcanc = 1.0;
-		goto blowup;
-		}
-*/
-	an += 1.0;
-	bn += 1.0;
-	n += 1.0;
-	}
-
-pdone:
-
-/* estimate error due to roundoff and cancellation */
-if( sum != 0.0 )
-	maxt /= fabsf(sum);
-maxt *= MACHEPF; 	/* this way avoids multiply overflow */
-pcanc = fabsf( MACHEPF * n  +  maxt );
-
-blowup:
-
-*err = pcanc;
-
-return( sum );
-}
-
-
-/*							hy1f1a()	*/
-/* asymptotic formula for hypergeometric function:
- *
- *        (    -a                         
- *  --    ( |z|                           
- * |  (b) ( -------- 2f0( a, 1+a-b, -1/x )
- *        (  --                           
- *        ( |  (b-a)                      
- *
- *
- *                                x    a-b                     )
- *                               e  |x|                        )
- *                             + -------- 2f0( b-a, 1-a, 1/x ) )
- *                                --                           )
- *                               |  (a)                        )
- */
-
-static float hy1f1af( float aa, float bb, float xx, float *err )
-{
-float a, b, x, h1, h2, t, u, temp, acanc, asum, err1, err2;
-
-a = aa;
-b = bb;
-x = xx;
-if( x == 0 )
-	{
-	acanc = 1.0;
-	asum = MAXNUMF;
-	goto adone;
-	}
-temp = logf( fabsf(x) );
-t = x + temp * (a-b);
-u = -temp * a;
-
-if( b > 0 )
-	{
-	temp = lgamf(b);
-	t += temp;
-	u += temp;
-	}
-
-h1 = hyp2f0f( a, a-b+1, -1.0/x, 1, &err1 );
-
-temp = expf(u) / gammaf(b-a);
-h1 *= temp;
-err1 *= temp;
-
-h2 = hyp2f0f( b-a, 1.0-a, 1.0/x, 2, &err2 );
-
-if( a < 0 )
-	temp = expf(t) / gammaf(a);
-else
-	temp = expf( t - lgamf(a) );
-
-h2 *= temp;
-err2 *= temp;
-
-if( x < 0.0 )
-	asum = h1;
-else
-	asum = h2;
-
-acanc = fabsf(err1) + fabsf(err2);
-
-
-if( b < 0 )
-	{
-	temp = gammaf(b);
-	asum *= temp;
-	acanc *= fabsf(temp);
-	}
-
-
-if( asum != 0.0 )
-	acanc /= fabsf(asum);
-
-acanc *= 30.0;	/* fudge factor, since error of asymptotic formula
-		 * often seems this much larger than advertised */
-
-adone:
-
-
-*err = acanc;
-return( asum );
-}
-
-/*							hyp2f0()	*/
-
-float hyp2f0f(float aa, float bb, float xx, int type, float *err)
-{
-float a, b, x, a0, alast, t, tlast, maxt;
-float n, an, bn, u, sum, temp;
-
-a = aa;
-b = bb;
-x = xx;
-an = a;
-bn = b;
-a0 = 1.0;
-alast = 1.0;
-sum = 0.0;
-n = 1.0;
-t = 1.0;
-tlast = 1.0e9;
-maxt = 0.0;
-
-do
-	{
-	if( an == 0 )
-		goto pdone;
-	if( bn == 0 )
-		goto pdone;
-
-	u = an * (bn * x / n);
-
-	/* check for blowup */
-	temp = fabsf(u);
-	if( (temp > 1.0 ) && (maxt > (MAXNUMF/temp)) )
-		goto error;
-
-	a0 *= u;
-	t = fabsf(a0);
-
-	/* terminating condition for asymptotic series */
-	if( t > tlast )
-		goto ndone;
-
-	tlast = t;
-	sum += alast;	/* the sum is one term behind */
-	alast = a0;
-
-	if( n > 200 )
-		goto ndone;
-
-	an += 1.0;
-	bn += 1.0;
-	n += 1.0;
-	if( t > maxt )
-		maxt = t;
-	}
-while( t > MACHEPF );
-
-
-pdone:	/* series converged! */
-
-/* estimate error due to roundoff and cancellation */
-*err = fabsf(  MACHEPF * (n + maxt)  );
-
-alast = a0;
-goto done;
-
-ndone:	/* series did not converge */
-
-/* The following "Converging factors" are supposed to improve accuracy,
- * but do not actually seem to accomplish very much. */
-
-n -= 1.0;
-x = 1.0/x;
-
-switch( type )	/* "type" given as subroutine argument */
-{
-case 1:
-	alast *= ( 0.5 + (0.125 + 0.25*b - 0.5*a + 0.25*x - 0.25*n)/x );
-	break;
-
-case 2:
-	alast *= 2.0/3.0 - b + 2.0*a + x - n;
-	break;
-
-default:
-	;
-}
-
-/* estimate error due to roundoff, cancellation, and nonconvergence */
-*err = MACHEPF * (n + maxt)  +  fabsf( a0 );
-
-
-done:
-sum += alast;
-return( sum );
-
-/* series blew up: */
-error:
-*err = MAXNUMF;
-mtherr( "hypergf", TLOSS );
-return( sum );
-}