summaryrefslogtreecommitdiff
path: root/libm/double/jv.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/double/jv.c')
-rw-r--r--libm/double/jv.c884
1 files changed, 884 insertions, 0 deletions
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 <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);
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-01 03:54:06 +0000