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

diff --git a/libm/float/struvef.c b/libm/float/struvef.c
deleted file mode 100644
index 4cf8854ed..000000000
--- a/libm/float/struvef.c
+++ /dev/null
@@ -1,315 +0,0 @@
-/*							struvef.c
- *
- *      Struve function
- *
- *
- *
- * SYNOPSIS:
- *
- * float v, x, y, struvef();
- *
- * y = struvef( v, x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Computes the Struve function Hv(x) of order v, argument x.
- * Negative x is rejected unless v is an integer.
- *
- * This module also contains the hypergeometric functions 1F2
- * and 3F0 and a routine for the Bessel function Yv(x) with
- * noninteger v.
- *
- *
- *
- * ACCURACY:
- *
- *  v varies from 0 to 10.
- *    Absolute error (relative error when |Hv(x)| > 1):
- * arithmetic   domain     # trials      peak         rms
- *    IEEE      -10,10      100000      9.0e-5      4.0e-6
- *
- */
-
-
-/*
-Cephes Math Library Release 2.2:  July, 1992
-Copyright 1984, 1987, 1989 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
-*/
-
-#include <math.h>
-#define DEBUG 0
-
-extern float MACHEPF, MAXNUMF, PIF;
-
-#define fabsf(x) ( (x) < 0 ? -(x) : (x) )
-
-#ifdef ANSIC
-float gammaf(float), powf(float, float), sqrtf(float);
-float yvf(float, float);
-float floorf(float), ynf(int, float);
-float jvf(float, float);
-float sinf(float), cosf(float);
-#else
-float gammaf(), powf(), sqrtf(), yvf();
-float floorf(), ynf(), jvf(), sinf(), cosf();
-#endif
-
-float onef2f( float aa, float bb, float cc, float xx, float *err )
-{
-float a, b, c, x, n, a0, sum, t;
-float an, bn, cn, max, z;
-
-a = aa;
-b = bb;
-c = cc;
-x = xx;
-an = a;
-bn = b;
-cn = c;
-a0 = 1.0;
-sum = 1.0;
-n = 1.0;
-t = 1.0;
-max = 0.0;
-
-do
-	{
-	if( an == 0 )
-		goto done;
-	if( bn == 0 )
-		goto error;
-	if( cn == 0 )
-		goto error;
-	if( (a0 > 1.0e34) || (n > 200) )
-		goto error;
-	a0 *= (an * x) / (bn * cn * n);
-	sum += a0;
-	an += 1.0;
-	bn += 1.0;
-	cn += 1.0;
-	n += 1.0;
-	z = fabsf( a0 );
-	if( z > max )
-		max = z;
-	if( sum != 0 )
-		t = fabsf( a0 / sum );
-	else
-		t = z;
-	}
-while( t > MACHEPF );
-
-done:
-
-*err = fabsf( MACHEPF*max /sum );
-
-#if DEBUG
-	printf(" onef2f cancellation error %.5E\n", *err );
-#endif
-
-goto xit;
-
-error:
-#if DEBUG
-printf("onef2f does not converge\n");
-#endif
-*err = MAXNUMF;
-
-xit:
-
-#if DEBUG
-printf("onef2( %.2E %.2E %.2E %.5E ) =  %.3E  %.6E\n", a, b, c, x, n, sum);
-#endif
-return(sum);
-}
-
-
-
-float threef0f( float aa, float bb, float cc, float xx, float *err )
-{
-float a, b, c, x, n, a0, sum, t, conv, conv1;
-float an, bn, cn, max, z;
-
-a = aa;
-b = bb;
-c = cc;
-x = xx;
-an = a;
-bn = b;
-cn = c;
-a0 = 1.0;
-sum = 1.0;
-n = 1.0;
-t = 1.0;
-max = 0.0;
-conv = 1.0e38;
-conv1 = conv;
-
-do
-	{
-	if( an == 0.0 )
-		goto done;
-	if( bn == 0.0 )
-		goto done;
-	if( cn == 0.0 )
-		goto done;
-	if( (a0 > 1.0e34) || (n > 200) )
-		goto error;
-	a0 *= (an * bn * cn * x) / n;
-	an += 1.0;
-	bn += 1.0;
-	cn += 1.0;
-	n += 1.0;
-	z = fabsf( a0 );
-	if( z > max )
-		max = z;
-	if( z >= conv )
-		{
-		if( (z < max) && (z > conv1) )
-			goto done;
-		}
-	conv1 = conv;
-	conv = z;
-	sum += a0;
-	if( sum != 0 )
-		t = fabsf( a0 / sum );
-	else
-		t = z;
-	}
-while( t > MACHEPF );
-
-done:
-
-t = fabsf( MACHEPF*max/sum );
-#if DEBUG
-	printf(" threef0f cancellation error %.5E\n", t );
-#endif
-
-max = fabsf( conv/sum );
-if( max > t )
-	t = max;
-#if DEBUG
-	printf(" threef0f convergence %.5E\n", max );
-#endif
-
-goto xit;
-
-error:
-#if DEBUG
-printf("threef0f does not converge\n");
-#endif
-t = MAXNUMF;
-
-xit:
-
-#if DEBUG
-printf("threef0f( %.2E %.2E %.2E %.5E ) =  %.3E  %.6E\n", a, b, c, x, n, sum);
-#endif
-
-*err = t;
-return(sum);
-}
-
-
-
-
-float struvef( float vv, float xx )
-{
-float v, x, y, ya, f, g, h, t;
-float onef2err, threef0err;
-
-v = vv;
-x = xx;
-f = floorf(v);
-if( (v < 0) && ( v-f == 0.5 ) )
-	{
-	y = jvf( -v, x );
-	f = 1.0 - f;
-	g =  2.0 * floorf(0.5*f);
-	if( g != f )
-		y = -y;
-	return(y);
-	}
-t = 0.25*x*x;
-f = fabsf(x);
-g = 1.5 * fabsf(v);
-if( (f > 30.0) && (f > g) )
-	{
-	onef2err = MAXNUMF;
-	y = 0.0;
-	}
-else
-	{
-	y = onef2f( 1.0, 1.5, 1.5+v, -t, &onef2err );
-	}
-
-if( (f < 18.0) || (x < 0.0) )
-	{
-	threef0err = MAXNUMF;
-	ya = 0.0;
-	}
-else
-	{
-	ya = threef0f( 1.0, 0.5, 0.5-v, -1.0/t, &threef0err );
-	}
-
-f = sqrtf( PIF );
-h = powf( 0.5*x, v-1.0 );
-
-if( onef2err <= threef0err )
-	{
-	g = gammaf( v + 1.5 );
-	y = y * h * t / ( 0.5 * f * g );
-	return(y);
-	}
-else
-	{
-	g = gammaf( v + 0.5 );
-	ya = ya * h / ( f * g );
-	ya = ya + yvf( v, x );
-	return(ya);
-	}
-}
-
-
-
-
-/* Bessel function of noninteger order
- */
-
-float yvf( float vv, float xx )
-{
-float v, x,  y, t;
-int n;
-
-v = vv;
-x = xx;
-y = floorf( v );
-if( y == v )
-	{
-	n = v;
-	y = ynf( n, x );
-	return( y );
-	}
-t = PIF * v;
-y = (cosf(t) * jvf( v, x ) - jvf( -v, x ))/sinf(t);
-return( y );
-}
-
-/* Crossover points between ascending series and asymptotic series
- * for Struve function
- *
- *	 v	 x
- * 
- *	 0	19.2
- *	 1	18.95
- *	 2	19.15
- *	 3	19.3
- *	 5	19.7
- *	10	21.35
- *	20	26.35
- *	30	32.31
- *	40	40.0
- */