uClibc now has a math library. muahahahaha!

 -Erik

1 files changed, 848 insertions, 0 deletions

diff --git a/libm/float/jvf.c b/libm/float/jvf.c
new file mode 100644
index 000000000..268a8e4eb
--- /dev/null
+++ b/libm/float/jvf.c
@@ -0,0 +1,848 @@
+/*							jvf.c
+ *
+ *	Bessel function of noninteger order
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float v, x, y, jvf();
+ *
+ * y = jvf( 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 single precision routine accepts negative v, but with
+ * reduced accuracy.
+ *
+ *
+ *
+ * ACCURACY:
+ * Results for integer v are indicated by *.
+ * Error criterion is absolute, except relative when |jv()| > 1.
+ *
+ * arithmetic     domain      # trials      peak         rms
+ *                v      x
+ *    IEEE       0,125  0,125   30000      2.0e-6      2.0e-7
+ *    IEEE     -17,0    0,125   30000      1.1e-5      4.0e-7
+ *    IEEE    -100,0    0,125    3000      1.5e-4      7.8e-6
+ */
+
+
+/*
+Cephes Math Library Release 2.2: June, 1992
+Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+
+#include <math.h>
+#define DEBUG 0
+
+extern float MAXNUMF, MACHEPF, MINLOGF, MAXLOGF, PIF;
+extern int sgngamf;
+
+/* BIG = 1/MACHEPF */
+#define BIG   16777216.
+
+#ifdef ANSIC
+float floorf(float), j0f(float), j1f(float);
+static float jnxf(float, float);
+static float jvsf(float, float);
+static float hankelf(float, float);
+static float jntf(float, float);
+static float recurf( float *, float, float * );
+float sqrtf(float), sinf(float), cosf(float);
+float lgamf(float), expf(float), logf(float), powf(float, float);
+float gammaf(float), cbrtf(float), acosf(float);
+int airyf(float, float *, float *, float *, float *);
+float polevlf(float, float *, int);
+#else
+float floorf(), j0f(), j1f();
+float sqrtf(), sinf(), cosf();
+float lgamf(), expf(), logf(), powf(), gammaf();
+float cbrtf(), polevlf(), acosf();
+void airyf();
+static float recurf(), jvsf(), hankelf(), jnxf(), jntf(), jvsf();
+#endif
+
+
+#define fabsf(x) ( (x) < 0 ? -(x) : (x) )
+
+float jvf( float nn, float xx )
+{
+float n, x, k, q, t, y, an, sign;
+int i, nint;
+
+n = nn;
+x = xx;
+nint = 0;	/* Flag for integer n */
+sign = 1.0;	/* Flag for sign inversion */
+an = fabsf( n );
+y = floorf( an );
+if( y == an )
+	{
+	nint = 1;
+	i = an - 16384.0 * floorf( 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( j0f(x) );
+	if( n == 1.0 )
+		return( sign * j1f(x) );
+	}
+
+if( (x < 0.0) && (y != an) )
+	{
+	mtherr( "jvf", DOMAIN );
+	y = 0.0;
+	goto done;
+	}
+
+y = fabsf(x);
+
+if( y < MACHEPF )
+	goto underf;
+
+/* Easy cases - x small compared to n */
+t = 3.6 * sqrtf(an);
+if( y < t )
+	return( sign * jvsf(n,x) );
+
+/* x large compared to n */
+k = 3.6 * sqrtf(y);
+if( (an < k) && (y > 6.0) )
+	return( sign * hankelf(n,x) );
+
+if( (n > -100) && (n < 14.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 = recurf( &n, x, &k );
+		if( k == 0.0 )
+			{
+			y = j0f(x)/q;
+			goto done;
+			}
+		if( k == 1.0 )
+			{
+			y = j1f(x)/q;
+			goto done;
+			}
+		}
+
+	if( n >= 0.0 )
+		{
+/* Recur backwards from a larger value of n
+ */
+		if( y > 1.3 * an )
+			goto recurdwn;
+		if( an > 1.3 * y )
+			goto recurdwn;
+		k = n;
+		y = 2.0*(y+an+1.0);
+		if( (y - n) > 33.0 )
+			y = n + 33.0;
+		y = n + floorf(y-n);
+		q = recurf( &y, x, &k );
+		y = jvsf(y,x) * q;
+		goto done;
+		}
+recurdwn:
+	if( an > (k + 3.0) )
+		{
+/* Recur backwards from n to k
+ */
+		if( n < 0.0 )
+			k = -k;
+		q = n - floorf(n);
+		k = floorf(k) + q;
+		if( n > 0.0 )
+			q = recurf( &n, x, &k );
+		else
+			{
+			t = k;
+			k = n;
+			q = recurf( &t, x, &k );
+			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
+ */
+ 	t = fabsf(k);
+	if( t < 26.0 )
+		t = (0.0083*t + 0.09)*t + 12.9;
+	else
+		t = 0.9 * t;
+
+	if( y > t ) /* y = |x| */
+		y = hankelf(k,x);
+	else
+		y = jvsf(k,x);
+#if DEBUG
+printf( "y = %.16e, q = %.16e\n", y, q );
+#endif
+	if( n > 0.0 )
+		y /= q;
+	else
+		y *= q;
+	}
+
+else
+	{
+/* For large positive 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( "jvf", TLOSS );
+		y = 0.0;
+		goto done;
+		}
+	t = y/an;
+	t /= an;
+	if( t > 0.3 )
+		y = hankelf(n,x);
+	else
+		y = jnxf(n,x);
+	}
+
+done:	return( sign * y);
+}
+
+/* Reduce the order by backward recurrence.
+ * AMS55 #9.1.27 and 9.1.73.
+ */
+
+static float recurf( float *n, float xx, float *newn )
+{
+float x, pkm2, pkm1, pk, pkp1, qkm2, qkm1;
+float k, ans, qk, xk, yk, r, t, kf, xinv;
+static float big = BIG;
+int nflag, ctr;
+
+x = xx;
+/* continued fraction for Jn(x)/Jn-1(x)  */
+if( *n < 0.0 )
+	nflag = 1;
+else
+	nflag = 0;
+
+fstart:
+
+#if DEBUG
+printf( "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 = fabsf( (ans - r)/r );
+		ans = r;
+		}
+	else
+		t = 1.0;
+
+	if( t < MACHEPF )
+		goto done;
+
+	if( fabsf(pk) > big )
+		{
+		pkm2 *= MACHEPF;
+		pkm1 *= MACHEPF;
+		qkm2 *= MACHEPF;
+		qkm1 *= MACHEPF;
+		}
+	}
+while( t > MACHEPF );
+
+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( fabsf(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;
+xinv = 1.0/x;
+do
+	{
+	pkm2 = (pkm1 * r  -  pk * x) * xinv;
+	pkp1 = pk;
+	pk = pkm1;
+	pkm1 = pkm2;
+	r -= 2.0;
+#if 0
+	t = fabsf(pkp1) + fabsf(pk);
+	if( (k > (kf + 2.5)) && (fabsf(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;
+		}
+#endif
+	k -= 1.0;
+	}
+while( k > (kf + 0.5) );
+
+#if 0
+/* Take the larger of the last two iterates
+ * on the theory that it may have less cancellation error.
+ */
+if( (kf >= 0.0) && (fabsf(pk) > fabsf(pkm1)) )
+	{
+	k += 1.0;
+	pkm2 = pk;
+	}
+#endif
+
+*newn = k;
+#if DEBUG
+printf( "newn %.6e\n", k );
+#endif
+return( pkm2 );
+}
+
+
+
+/* Ascending power series for Jv(x).
+ * AMS55 #9.1.10.
+ */
+
+static float jvsf( float nn, float xx )
+{
+float n, x, t, u, y, z, k, ay;
+
+#if DEBUG
+printf( "jvsf: " );
+#endif
+n = nn;
+x = xx;
+z = -0.25 * x * x;
+u = 1.0;
+y = u;
+k = 1.0;
+t = 1.0;
+
+while( t > MACHEPF )
+	{
+	u *= z / (k * (n+k));
+	y += u;
+	k += 1.0;
+	t = fabsf(u);
+	if( (ay = fabsf(y)) > 1.0 )
+		t /= ay;
+	}
+
+if( x < 0.0 )
+	{
+	y = y * powf( 0.5 * x, n ) / gammaf( n + 1.0 );
+	}
+else
+	{
+	t = n * logf(0.5*x) - lgamf(n + 1.0);
+	if( t < -MAXLOGF )
+		{
+		return( 0.0 );
+		}
+	if( t > MAXLOGF )
+		{
+		t = logf(y) + t;
+		if( t > MAXLOGF )
+			{
+			mtherr( "jvf", OVERFLOW );
+			return( MAXNUMF );
+			}
+		else
+			{
+			y = sgngamf * expf(t);
+			return(y);
+			}
+		}
+	y = sgngamf * y * expf( t );
+	}
+#if DEBUG
+printf( "y = %.8e\n", y );
+#endif
+return(y);
+}
+
+/* Hankel's asymptotic expansion
+ * for large x.
+ * AMS55 #9.2.5.
+ */
+static float hankelf( float nn, float xx )
+{
+float n, x, t, u, z, k, sign, conv;
+float p, q, j, m, pp, qq;
+int flag;
+
+#if DEBUG
+printf( "hankelf: " );
+#endif
+n = nn;
+x = xx;
+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 > MACHEPF )
+	{
+	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 = fabsf(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) * PIF;
+t = sqrtf( 2.0/(PIF*x) ) * ( pp * cosf(u) - qq * sinf(u) );
+return( t );
+}
+
+
+/* Asymptotic expansion for large n.
+ * AMS55 #9.3.35.
+ */
+
+static float 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 float 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 float P1[] = {
+ -2.083333333333333333333333E-1,
+  1.250000000000000000000000E-1
+};
+static float P2[] = {
+  3.342013888888888888888889E-1,
+ -4.010416666666666666666667E-1,
+  7.031250000000000000000000E-2
+};
+static float P3[] = {
+ -1.025812596450617283950617E+0,
+  1.846462673611111111111111E+0,
+ -8.912109375000000000000000E-1,
+  7.324218750000000000000000E-2
+};
+static float P4[] = {
+  4.669584423426247427983539E+0,
+ -1.120700261622299382716049E+1,
+  8.789123535156250000000000E+0,
+ -2.364086914062500000000000E+0,
+  1.121520996093750000000000E-1
+};
+static float P5[] = {
+ -2.8212072558200244877E1,
+  8.4636217674600734632E1,
+ -9.1818241543240017361E1,
+  4.2534998745388454861E1,
+ -7.3687943594796316964E0,
+  2.27108001708984375E-1
+};
+static float P6[] = {
+  2.1257013003921712286E2,
+ -7.6525246814118164230E2,
+  1.0599904525279998779E3,
+ -6.9957962737613254123E2,
+  2.1819051174421159048E2,
+ -2.6491430486951555525E1,
+  5.7250142097473144531E-1
+};
+static float P7[] = {
+ -1.9194576623184069963E3,
+  8.0617221817373093845E3,
+ -1.3586550006434137439E4,
+  1.1655393336864533248E4,
+ -5.3056469786134031084E3,
+  1.2009029132163524628E3,
+ -1.0809091978839465550E2,
+  1.7277275025844573975E0
+};
+
+
+static float jnxf( float nn, float xx )
+{
+float n, x, zeta, sqz, zz, zp, np;
+float cbn, n23, t, z, sz;
+float pp, qq, z32i, zzi;
+float ak, bk, akl, bkl;
+int sign, doa, dob, nflg, k, s, tk, tkp1, m;
+static float u[8];
+static float ai, aip, bi, bip;
+
+n = nn;
+x = xx;
+/* Test for x very close to n.
+ * Use expansion for transition region if so.
+ */
+cbn = cbrtf(n);
+z = (x - n)/cbn;
+if( (fabsf(z) <= 0.7) || (n < 0.0) )
+	return( jntf(n,x) );
+z = x/n;
+zz = 1.0 - z*z;
+if( zz == 0.0 )
+	return(0.0);
+
+if( zz > 0.0 )
+	{
+	sz = sqrtf( zz );
+	t = 1.5 * (logf( (1.0+sz)/z ) - sz );	/* zeta ** 3/2		*/
+	zeta = cbrtf( t * t );
+	nflg = 1;
+	}
+else
+	{
+	sz = sqrtf(-zz);
+	t = 1.5 * (sz - acosf(1.0/z));
+	zeta = -cbrtf( t * t );
+	nflg = -1;
+	}
+z32i = fabsf(1.0/t);
+sqz = cbrtf(t);
+
+/* Airy function */
+n23 = cbrtf( n * n );
+t = n23 * zeta;
+
+#if DEBUG
+printf("zeta %.5E, Airyf(%.5E)\n", zeta, t );
+#endif
+airyf( t, &ai, &aip, &bi, &bip );
+
+/* polynomials in expansion */
+u[0] = 1.0;
+zzi = 1.0/zz;
+u[1] = polevlf( zzi, P1, 1 )/sz;
+u[2] = polevlf( zzi, P2, 2 )/zz;
+u[3] = polevlf( zzi, P3, 3 )/(sz*zz);
+pp = zz*zz;
+u[4] = polevlf( zzi, P4, 4 )/pp;
+u[5] = polevlf( zzi, P5, 5 )/(pp*sz);
+pp *= zz;
+u[6] = polevlf( zzi, P6, 6 )/pp;
+u[7] = polevlf( 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 = MAXNUMF;
+bkl = MAXNUMF;
+
+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 = fabsf(ak);
+		if( t < akl )
+			{
+			akl = t;
+			pp += ak;
+			}
+		else
+			doa = 0;
+		}
+
+	if( dob )
+		{
+		bk += lambda[tkp1] * zp * u[0];
+		bk *= -np/sqz;
+		t = fabsf(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 < MACHEPF )
+		break;
+	np /= n*n;
+	}
+
+/* normalizing factor ( 4*zeta/(1 - z**2) )**1/4	*/
+t = 4.0 * zeta/zz;
+t = sqrtf( sqrtf(t) );
+
+t *= ai*pp/cbrtf(n)  +  aip*qq/(n23*n);
+return(t);
+}
+
+/* Asymptotic expansion for transition region,
+ * n large and x close to n.
+ * AMS55 #9.3.23.
+ */
+
+static float PF2[] = {
+ -9.0000000000000000000e-2,
+  8.5714285714285714286e-2
+};
+static float PF3[] = {
+  1.3671428571428571429e-1,
+ -5.4920634920634920635e-2,
+ -4.4444444444444444444e-3
+};
+static float PF4[] = {
+  1.3500000000000000000e-3,
+ -1.6036054421768707483e-1,
+  4.2590187590187590188e-2,
+  2.7330447330447330447e-3
+};
+static float PG1[] = {
+ -2.4285714285714285714e-1,
+  1.4285714285714285714e-2
+};
+static float PG2[] = {
+ -9.0000000000000000000e-3,
+  1.9396825396825396825e-1,
+ -1.1746031746031746032e-2
+};
+static float PG3[] = {
+  1.9607142857142857143e-2,
+ -1.5983694083694083694e-1,
+  6.3838383838383838384e-3
+};
+
+
+static float jntf( float nn, float xx )
+{
+float n, x, z, zz, z3;
+float cbn, n23, cbtwo;
+float ai, aip, bi, bip;	/* Airy functions */
+float nk, fk, gk, pp, qq;
+float F[5], G[4];
+int k;
+
+n = nn;
+x = xx;
+cbn = cbrtf(n);
+z = (x - n)/cbn;
+cbtwo = cbrtf( 2.0 );
+
+/* Airy function */
+zz = -cbtwo * z;
+airyf( 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] = polevlf( z3, PF2, 1 ) * zz;
+F[3] = polevlf( z3, PF3, 2 );
+F[4] = polevlf( z3, PF4, 3 ) * z;
+G[0] = 0.3 * zz;
+G[1] = polevlf( z3, PG1, 1 );
+G[2] = polevlf( z3, PG2, 2 ) * z;
+G[3] = polevlf( 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 = cbrtf( 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  +  cbrtf(4.0) * aip * qq/n;
+return(fk);
+}