From 1077fa4d772832f77a677ce7fb7c2d513b959e3f Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 10 May 2001 00:40:28 +0000 Subject: uClibc now has a math library. muahahahaha! -Erik --- libm/double/jv.c | 884 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 884 insertions(+) create mode 100644 libm/double/jv.c (limited to 'libm/double/jv.c') diff --git a/libm/double/jv.c b/libm/double/jv.c new file mode 100644 index 000000000..5b8af3663 --- /dev/null +++ b/libm/double/jv.c @@ -0,0 +1,884 @@ +/* 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 +#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); +} -- cgit v1.2.3