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

diff --git a/libm/double/jv.c b/libm/double/jv.c
deleted file mode 100644
index 5b8af3663..000000000
--- a/libm/double/jv.c
+++ /dev/null
@@ -1,884 +0,0 @@
-/*							jv.c
- *
- *	Bessel function of noninteger order
- *
- *
- *
- * SYNOPSIS:
- *
- * double v, x, y, jv();
- *
- * y = jv( v, x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns Bessel function of order v of the argument,
- * where v is real.  Negative x is allowed if v is an integer.
- *
- * Several expansions are included: the ascending power
- * series, the Hankel expansion, and two transitional
- * expansions for large v.  If v is not too large, it
- * is reduced by recurrence to a region of best accuracy.
- * The transitional expansions give 12D accuracy for v > 500.
- *
- *
- *
- * ACCURACY:
- * Results for integer v are indicated by *, where x and v
- * both vary from -125 to +125.  Otherwise,
- * x ranges from 0 to 125, v ranges as indicated by "domain."
- * Error criterion is absolute, except relative when |jv()| > 1.
- *
- * arithmetic  v domain  x domain    # trials      peak       rms
- *    IEEE      0,125     0,125      100000      4.6e-15    2.2e-16
- *    IEEE   -125,0       0,125       40000      5.4e-11    3.7e-13
- *    IEEE      0,500     0,500       20000      4.4e-15    4.0e-16
- * Integer v:
- *    IEEE   -125,125   -125,125      50000      3.5e-15*   1.9e-16*
- *
- */
-
-
-/*
-Cephes Math Library Release 2.8:  June, 2000
-Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
-*/
-
-
-#include <math.h>
-#define DEBUG 0
-
-#ifdef DEC
-#define MAXGAM 34.84425627277176174
-#else
-#define MAXGAM 171.624376956302725
-#endif
-
-#ifdef ANSIPROT
-extern int airy ( double, double *, double *, double *, double * );
-extern double fabs ( double );
-extern double floor ( double );
-extern double frexp ( double, int * );
-extern double polevl ( double, void *, int );
-extern double j0 ( double );
-extern double j1 ( double );
-extern double sqrt ( double );
-extern double cbrt ( double );
-extern double exp ( double );
-extern double log ( double );
-extern double sin ( double );
-extern double cos ( double );
-extern double acos ( double );
-extern double pow ( double, double );
-extern double gamma ( double );
-extern double lgam ( double );
-static double recur(double *, double, double *, int);
-static double jvs(double, double);
-static double hankel(double, double);
-static double jnx(double, double);
-static double jnt(double, double);
-#else
-int airy();
-double fabs(), floor(), frexp(), polevl(), j0(), j1(), sqrt(), cbrt();
-double exp(), log(), sin(), cos(), acos(), pow(), gamma(), lgam();
-static double recur(), jvs(), hankel(), jnx(), jnt();
-#endif
-
-extern double MAXNUM, MACHEP, MINLOG, MAXLOG;
-#define BIG  1.44115188075855872E+17
-
-double jv( n, x )
-double n, x;
-{
-double k, q, t, y, an;
-int i, sign, nint;
-
-nint = 0;	/* Flag for integer n */
-sign = 1;	/* Flag for sign inversion */
-an = fabs( n );
-y = floor( an );
-if( y == an )
-	{
-	nint = 1;
-	i = an - 16384.0 * floor( an/16384.0 );
-	if( n < 0.0 )
-		{
-		if( i & 1 )
-			sign = -sign;
-		n = an;
-		}
-	if( x < 0.0 )
-		{
-		if( i & 1 )
-			sign = -sign;
-		x = -x;
-		}
-	if( n == 0.0 )
-		return( j0(x) );
-	if( n == 1.0 )
-		return( sign * j1(x) );
-	}
-
-if( (x < 0.0) && (y != an) )
-	{
-	mtherr( "Jv", DOMAIN );
-	y = 0.0;
-	goto done;
- 	}
-
-y = fabs(x);
-
-if( y < MACHEP )
-	goto underf;
-
-k = 3.6 * sqrt(y);
-t = 3.6 * sqrt(an);
-if( (y < t) && (an > 21.0) )
-	return( sign * jvs(n,x) );
-if( (an < k) && (y > 21.0) )
-	return( sign * hankel(n,x) );
-
-if( an < 500.0 )
-	{
-/* Note: if x is too large, the continued
- * fraction will fail; but then the
- * Hankel expansion can be used.
- */
-	if( nint != 0 )
-		{
-		k = 0.0;
-		q = recur( &n, x, &k, 1 );
-		if( k == 0.0 )
-			{
-			y = j0(x)/q;
-			goto done;
-			}
-		if( k == 1.0 )
-			{
-			y = j1(x)/q;
-			goto done;
-			}
-		}
-
-if( an > 2.0 * y )
-	goto rlarger;
-
-	if( (n >= 0.0) && (n < 20.0)
-		&& (y > 6.0) && (y < 20.0) )
-		{
-/* Recur backwards from a larger value of n
- */
-rlarger:
-		k = n;
-
-		y = y + an + 1.0;
-		if( y < 30.0 )
-			y = 30.0;
-		y = n + floor(y-n);
-		q = recur( &y, x, &k, 0 );
-		y = jvs(y,x) * q;
-		goto done;
-		}
-
-	if( k <= 30.0 )
-		{
-		k = 2.0;
-		}
-	else if( k < 90.0 )
-		{
-		k = (3*k)/4;
-		}
-	if( an > (k + 3.0) )
-		{
-		if( n < 0.0 )
-			k = -k;
-		q = n - floor(n);
-		k = floor(k) + q;
-		if( n > 0.0 )
-			q = recur( &n, x, &k, 1 );
-		else
-			{
-			t = k;
-			k = n;
-			q = recur( &t, x, &k, 1 );
-			k = t;
-			}
-		if( q == 0.0 )
-			{
-underf:
-			y = 0.0;
-			goto done;
-			}
-		}
-	else
-		{
-		k = n;
-		q = 1.0;
-		}
-
-/* boundary between convergence of
- * power series and Hankel expansion
- */
-	y = fabs(k);
-	if( y < 26.0 )
-		t = (0.0083*y + 0.09)*y + 12.9;
-	else
-		t = 0.9 * y;
-
-	if( x > t )
-		y = hankel(k,x);
-	else
-		y = jvs(k,x);
-#if DEBUG
-printf( "y = %.16e, recur q = %.16e\n", y, q );
-#endif
-	if( n > 0.0 )
-		y /= q;
-	else
-		y *= q;
-	}
-
-else
-	{
-/* For large n, use the uniform expansion
- * or the transitional expansion.
- * But if x is of the order of n**2,
- * these may blow up, whereas the
- * Hankel expansion will then work.
- */
-	if( n < 0.0 )
-		{
-		mtherr( "Jv", TLOSS );
-		y = 0.0;
-		goto done;
-		}
-	t = x/n;
-	t /= n;
-	if( t > 0.3 )
-		y = hankel(n,x);
-	else
-		y = jnx(n,x);
-	}
-
-done:	return( sign * y);
-}
-
-/* Reduce the order by backward recurrence.
- * AMS55 #9.1.27 and 9.1.73.
- */
-
-static double recur( n, x, newn, cancel )
-double *n;
-double x;
-double *newn;
-int cancel;
-{
-double pkm2, pkm1, pk, qkm2, qkm1;
-/* double pkp1; */
-double k, ans, qk, xk, yk, r, t, kf;
-static double big = BIG;
-int nflag, ctr;
-
-/* continued fraction for Jn(x)/Jn-1(x)  */
-if( *n < 0.0 )
-	nflag = 1;
-else
-	nflag = 0;
-
-fstart:
-
-#if DEBUG
-printf( "recur: n = %.6e, newn = %.6e, cfrac = ", *n, *newn );
-#endif
-
-pkm2 = 0.0;
-qkm2 = 1.0;
-pkm1 = x;
-qkm1 = *n + *n;
-xk = -x * x;
-yk = qkm1;
-ans = 1.0;
-ctr = 0;
-do
-	{
-	yk += 2.0;
-	pk = pkm1 * yk +  pkm2 * xk;
-	qk = qkm1 * yk +  qkm2 * xk;
-	pkm2 = pkm1;
-	pkm1 = pk;
-	qkm2 = qkm1;
-	qkm1 = qk;
-	if( qk != 0 )
-		r = pk/qk;
-	else
-		r = 0.0;
-	if( r != 0 )
-		{
-		t = fabs( (ans - r)/r );
-		ans = r;
-		}
-	else
-		t = 1.0;
-
-	if( ++ctr > 1000 )
-		{
-		mtherr( "jv", UNDERFLOW );
-		goto done;
-		}
-	if( t < MACHEP )
-		goto done;
-
-	if( fabs(pk) > big )
-		{
-		pkm2 /= big;
-		pkm1 /= big;
-		qkm2 /= big;
-		qkm1 /= big;
-		}
-	}
-while( t > MACHEP );
-
-done:
-
-#if DEBUG
-printf( "%.6e\n", ans );
-#endif
-
-/* Change n to n-1 if n < 0 and the continued fraction is small
- */
-if( nflag > 0 )
-	{
-	if( fabs(ans) < 0.125 )
-		{
-		nflag = -1;
-		*n = *n - 1.0;
-		goto fstart;
-		}
-	}
-
-
-kf = *newn;
-
-/* backward recurrence
- *              2k
- *  J   (x)  =  --- J (x)  -  J   (x)
- *   k-1         x   k         k+1
- */
-
-pk = 1.0;
-pkm1 = 1.0/ans;
-k = *n - 1.0;
-r = 2 * k;
-do
-	{
-	pkm2 = (pkm1 * r  -  pk * x) / x;
-	/*	pkp1 = pk; */
-	pk = pkm1;
-	pkm1 = pkm2;
-	r -= 2.0;
-/*
-	t = fabs(pkp1) + fabs(pk);
-	if( (k > (kf + 2.5)) && (fabs(pkm1) < 0.25*t) )
-		{
-		k -= 1.0;
-		t = x*x;
-		pkm2 = ( (r*(r+2.0)-t)*pk - r*x*pkp1 )/t;
-		pkp1 = pk;
-		pk = pkm1;
-		pkm1 = pkm2;
-		r -= 2.0;
-		}
-*/
-	k -= 1.0;
-	}
-while( k > (kf + 0.5) );
-
-/* Take the larger of the last two iterates
- * on the theory that it may have less cancellation error.
- */
-
-if( cancel )
-	{
-	if( (kf >= 0.0) && (fabs(pk) > fabs(pkm1)) )
-		{
-		k += 1.0;
-		pkm2 = pk;
-		}
-	}
-*newn = k;
-#if DEBUG
-printf( "newn %.6e rans %.6e\n", k, pkm2 );
-#endif
-return( pkm2 );
-}
-
-
-
-/* Ascending power series for Jv(x).
- * AMS55 #9.1.10.
- */
-
-extern double PI;
-extern int sgngam;
-
-static double jvs( n, x )
-double n, x;
-{
-double t, u, y, z, k;
-int ex;
-
-z = -x * x / 4.0;
-u = 1.0;
-y = u;
-k = 1.0;
-t = 1.0;
-
-while( t > MACHEP )
-	{
-	u *= z / (k * (n+k));
-	y += u;
-	k += 1.0;
-	if( y != 0 )
-		t = fabs( u/y );
-	}
-#if DEBUG
-printf( "power series=%.5e ", y );
-#endif
-t = frexp( 0.5*x, &ex );
-ex = ex * n;
-if(  (ex > -1023)
-  && (ex < 1023) 
-  && (n > 0.0)
-  && (n < (MAXGAM-1.0)) )
-	{
-	t = pow( 0.5*x, n ) / gamma( n + 1.0 );
-#if DEBUG
-printf( "pow(.5*x, %.4e)/gamma(n+1)=%.5e\n", n, t );
-#endif
-	y *= t;
-	}
-else
-	{
-#if DEBUG
-	z = n * log(0.5*x);
-	k = lgam( n+1.0 );
-	t = z - k;
-	printf( "log pow=%.5e, lgam(%.4e)=%.5e\n", z, n+1.0, k );
-#else
-	t = n * log(0.5*x) - lgam(n + 1.0);
-#endif
-	if( y < 0 )
-		{
-		sgngam = -sgngam;
-		y = -y;
-		}
-	t += log(y);
-#if DEBUG
-printf( "log y=%.5e\n", log(y) );
-#endif
-	if( t < -MAXLOG )
-		{
-		return( 0.0 );
-		}
-	if( t > MAXLOG )
-		{
-		mtherr( "Jv", OVERFLOW );
-		return( MAXNUM );
-		}
-	y = sgngam * exp( t );
-	}
-return(y);
-}
-
-/* Hankel's asymptotic expansion
- * for large x.
- * AMS55 #9.2.5.
- */
-
-static double hankel( n, x )
-double n, x;
-{
-double t, u, z, k, sign, conv;
-double p, q, j, m, pp, qq;
-int flag;
-
-m = 4.0*n*n;
-j = 1.0;
-z = 8.0 * x;
-k = 1.0;
-p = 1.0;
-u = (m - 1.0)/z;
-q = u;
-sign = 1.0;
-conv = 1.0;
-flag = 0;
-t = 1.0;
-pp = 1.0e38;
-qq = 1.0e38;
-
-while( t > MACHEP )
-	{
-	k += 2.0;
-	j += 1.0;
-	sign = -sign;
-	u *= (m - k * k)/(j * z);
-	p += sign * u;
-	k += 2.0;
-	j += 1.0;
-	u *= (m - k * k)/(j * z);
-	q += sign * u;
-	t = fabs(u/p);
-	if( t < conv )
-		{
-		conv = t;
-		qq = q;
-		pp = p;
-		flag = 1;
-		}
-/* stop if the terms start getting larger */
-	if( (flag != 0) && (t > conv) )
-		{
-#if DEBUG
-		printf( "Hankel: convergence to %.4E\n", conv );
-#endif
-		goto hank1;
-		}
-	}	
-
-hank1:
-u = x - (0.5*n + 0.25) * PI;
-t = sqrt( 2.0/(PI*x) ) * ( pp * cos(u) - qq * sin(u) );
-#if DEBUG
-printf( "hank: %.6e\n", t );
-#endif
-return( t );
-}
-
-
-/* Asymptotic expansion for large n.
- * AMS55 #9.3.35.
- */
-
-static double lambda[] = {
-  1.0,
-  1.041666666666666666666667E-1,
-  8.355034722222222222222222E-2,
-  1.282265745563271604938272E-1,
-  2.918490264641404642489712E-1,
-  8.816272674437576524187671E-1,
-  3.321408281862767544702647E+0,
-  1.499576298686255465867237E+1,
-  7.892301301158651813848139E+1,
-  4.744515388682643231611949E+2,
-  3.207490090890661934704328E+3
-};
-static double mu[] = {
-  1.0,
- -1.458333333333333333333333E-1,
- -9.874131944444444444444444E-2,
- -1.433120539158950617283951E-1,
- -3.172272026784135480967078E-1,
- -9.424291479571202491373028E-1,
- -3.511203040826354261542798E+0,
- -1.572726362036804512982712E+1,
- -8.228143909718594444224656E+1,
- -4.923553705236705240352022E+2,
- -3.316218568547972508762102E+3
-};
-static double P1[] = {
- -2.083333333333333333333333E-1,
-  1.250000000000000000000000E-1
-};
-static double P2[] = {
-  3.342013888888888888888889E-1,
- -4.010416666666666666666667E-1,
-  7.031250000000000000000000E-2
-};
-static double P3[] = {
- -1.025812596450617283950617E+0,
-  1.846462673611111111111111E+0,
- -8.912109375000000000000000E-1,
-  7.324218750000000000000000E-2
-};
-static double P4[] = {
-  4.669584423426247427983539E+0,
- -1.120700261622299382716049E+1,
-  8.789123535156250000000000E+0,
- -2.364086914062500000000000E+0,
-  1.121520996093750000000000E-1
-};
-static double P5[] = {
- -2.8212072558200244877E1,
-  8.4636217674600734632E1,
- -9.1818241543240017361E1,
-  4.2534998745388454861E1,
- -7.3687943594796316964E0,
-  2.27108001708984375E-1
-};
-static double P6[] = {
-  2.1257013003921712286E2,
- -7.6525246814118164230E2,
-  1.0599904525279998779E3,
- -6.9957962737613254123E2,
-  2.1819051174421159048E2,
- -2.6491430486951555525E1,
-  5.7250142097473144531E-1
-};
-static double P7[] = {
- -1.9194576623184069963E3,
-  8.0617221817373093845E3,
- -1.3586550006434137439E4,
-  1.1655393336864533248E4,
- -5.3056469786134031084E3,
-  1.2009029132163524628E3,
- -1.0809091978839465550E2,
-  1.7277275025844573975E0
-};
-
-
-static double jnx( n, x )
-double n, x;
-{
-double zeta, sqz, zz, zp, np;
-double cbn, n23, t, z, sz;
-double pp, qq, z32i, zzi;
-double ak, bk, akl, bkl;
-int sign, doa, dob, nflg, k, s, tk, tkp1, m;
-static double u[8];
-static double ai, aip, bi, bip;
-
-/* Test for x very close to n.
- * Use expansion for transition region if so.
- */
-cbn = cbrt(n);
-z = (x - n)/cbn;
-if( fabs(z) <= 0.7 )
-	return( jnt(n,x) );
-
-z = x/n;
-zz = 1.0 - z*z;
-if( zz == 0.0 )
-	return(0.0);
-
-if( zz > 0.0 )
-	{
-	sz = sqrt( zz );
-	t = 1.5 * (log( (1.0+sz)/z ) - sz );	/* zeta ** 3/2		*/
-	zeta = cbrt( t * t );
-	nflg = 1;
-	}
-else
-	{
-	sz = sqrt(-zz);
-	t = 1.5 * (sz - acos(1.0/z));
-	zeta = -cbrt( t * t );
-	nflg = -1;
-	}
-z32i = fabs(1.0/t);
-sqz = cbrt(t);
-
-/* Airy function */
-n23 = cbrt( n * n );
-t = n23 * zeta;
-
-#if DEBUG
-printf("zeta %.5E, Airy(%.5E)\n", zeta, t );
-#endif
-airy( t, &ai, &aip, &bi, &bip );
-
-/* polynomials in expansion */
-u[0] = 1.0;
-zzi = 1.0/zz;
-u[1] = polevl( zzi, P1, 1 )/sz;
-u[2] = polevl( zzi, P2, 2 )/zz;
-u[3] = polevl( zzi, P3, 3 )/(sz*zz);
-pp = zz*zz;
-u[4] = polevl( zzi, P4, 4 )/pp;
-u[5] = polevl( zzi, P5, 5 )/(pp*sz);
-pp *= zz;
-u[6] = polevl( zzi, P6, 6 )/pp;
-u[7] = polevl( zzi, P7, 7 )/(pp*sz);
-
-#if DEBUG
-for( k=0; k<=7; k++ )
-	printf( "u[%d] = %.5E\n", k, u[k] );
-#endif
-
-pp = 0.0;
-qq = 0.0;
-np = 1.0;
-/* flags to stop when terms get larger */
-doa = 1;
-dob = 1;
-akl = MAXNUM;
-bkl = MAXNUM;
-
-for( k=0; k<=3; k++ )
-	{
-	tk = 2 * k;
-	tkp1 = tk + 1;
-	zp = 1.0;
-	ak = 0.0;
-	bk = 0.0;
-	for( s=0; s<=tk; s++ )
-		{
-		if( doa )
-			{
-			if( (s & 3) > 1 )
-				sign = nflg;
-			else
-				sign = 1;
-			ak += sign * mu[s] * zp * u[tk-s];
-			}
-
-		if( dob )
-			{
-			m = tkp1 - s;
-			if( ((m+1) & 3) > 1 )
-				sign = nflg;
-			else
-				sign = 1;
-			bk += sign * lambda[s] * zp * u[m];
-			}
-		zp *= z32i;
-		}
-
-	if( doa )
-		{
-		ak *= np;
-		t = fabs(ak);
-		if( t < akl )
-			{
-			akl = t;
-			pp += ak;
-			}
-		else
-			doa = 0;
-		}
-
-	if( dob )
-		{
-		bk += lambda[tkp1] * zp * u[0];
-		bk *= -np/sqz;
-		t = fabs(bk);
-		if( t < bkl )
-			{
-			bkl = t;
-			qq += bk;
-			}
-		else
-			dob = 0;
-		}
-#if DEBUG
-	printf("a[%d] %.5E, b[%d] %.5E\n", k, ak, k, bk );
-#endif
-	if( np < MACHEP )
-		break;
-	np /= n*n;
-	}
-
-/* normalizing factor ( 4*zeta/(1 - z**2) )**1/4	*/
-t = 4.0 * zeta/zz;
-t = sqrt( sqrt(t) );
-
-t *= ai*pp/cbrt(n)  +  aip*qq/(n23*n);
-return(t);
-}
-
-/* Asymptotic expansion for transition region,
- * n large and x close to n.
- * AMS55 #9.3.23.
- */
-
-static double PF2[] = {
- -9.0000000000000000000e-2,
-  8.5714285714285714286e-2
-};
-static double PF3[] = {
-  1.3671428571428571429e-1,
- -5.4920634920634920635e-2,
- -4.4444444444444444444e-3
-};
-static double PF4[] = {
-  1.3500000000000000000e-3,
- -1.6036054421768707483e-1,
-  4.2590187590187590188e-2,
-  2.7330447330447330447e-3
-};
-static double PG1[] = {
- -2.4285714285714285714e-1,
-  1.4285714285714285714e-2
-};
-static double PG2[] = {
- -9.0000000000000000000e-3,
-  1.9396825396825396825e-1,
- -1.1746031746031746032e-2
-};
-static double PG3[] = {
-  1.9607142857142857143e-2,
- -1.5983694083694083694e-1,
-  6.3838383838383838384e-3
-};
-
-
-static double jnt( n, x )
-double n, x;
-{
-double z, zz, z3;
-double cbn, n23, cbtwo;
-double ai, aip, bi, bip;	/* Airy functions */
-double nk, fk, gk, pp, qq;
-double F[5], G[4];
-int k;
-
-cbn = cbrt(n);
-z = (x - n)/cbn;
-cbtwo = cbrt( 2.0 );
-
-/* Airy function */
-zz = -cbtwo * z;
-airy( zz, &ai, &aip, &bi, &bip );
-
-/* polynomials in expansion */
-zz = z * z;
-z3 = zz * z;
-F[0] = 1.0;
-F[1] = -z/5.0;
-F[2] = polevl( z3, PF2, 1 ) * zz;
-F[3] = polevl( z3, PF3, 2 );
-F[4] = polevl( z3, PF4, 3 ) * z;
-G[0] = 0.3 * zz;
-G[1] = polevl( z3, PG1, 1 );
-G[2] = polevl( z3, PG2, 2 ) * z;
-G[3] = polevl( z3, PG3, 2 ) * zz;
-#if DEBUG
-for( k=0; k<=4; k++ )
-	printf( "F[%d] = %.5E\n", k, F[k] );
-for( k=0; k<=3; k++ )
-	printf( "G[%d] = %.5E\n", k, G[k] );
-#endif
-pp = 0.0;
-qq = 0.0;
-nk = 1.0;
-n23 = cbrt( n * n );
-
-for( k=0; k<=4; k++ )
-	{
-	fk = F[k]*nk;
-	pp += fk;
-	if( k != 4 )
-		{
-		gk = G[k]*nk;
-		qq += gk;
-		}
-#if DEBUG
-	printf("fk[%d] %.5E, gk[%d] %.5E\n", k, fk, k, gk );
-#endif
-	nk /= n23;
-	}
-
-fk = cbtwo * ai * pp/cbn  +  cbrt(4.0) * aip * qq/n;
-return(fk);
-}