lib/libm/jvf.c

   1 /*                                                      jvf.c
   2  *
   3  *      Bessel function of noninteger order
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float v, x, y, jvf();
  10  *
  11  * y = jvf( v, x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns Bessel function of order v of the argument,
  18  * where v is real.  Negative x is allowed if v is an integer.
  19  *
  20  * Several expansions are included: the ascending power
  21  * series, the Hankel expansion, and two transitional
  22  * expansions for large v.  If v is not too large, it
  23  * is reduced by recurrence to a region of best accuracy.
  24  *
  25  * The single precision routine accepts negative v, but with
  26  * reduced accuracy.
  27  *
  28  *
  29  *
  30  * ACCURACY:
  31  * Results for integer v are indicated by *.
  32  * Error criterion is absolute, except relative when |jv()| > 1.
  33  *
  34  * arithmetic     domain      # trials      peak         rms
  35  *                v      x
  36  *    IEEE       0,125  0,125   30000      2.0e-6      2.0e-7
  37  *    IEEE     -17,0    0,125   30000      1.1e-5      4.0e-7
  38  *    IEEE    -100,0    0,125    3000      1.5e-4      7.8e-6
  39  */
  40 \f
  41
  42 /*
  43 Cephes Math Library Release 2.2: June, 1992
  44 Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
  45 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  46 */
  47
  48
  49 #include "mconf.h"
  50 #define DEBUG 0
  51
  52 extern float MAXNUMF, MACHEPF, MINLOGF, MAXLOGF, PIF;
  53 extern int sgngamf;
  54
  55 /* BIG = 1/MACHEPF */
  56 #define BIG   16777216.
  57
  58 #ifdef ANSIC
  59 float floorf(float), j0f(float), j1f(float);
  60 static float jnxf(float, float);
  61 static float jvsf(float, float);
  62 static float hankelf(float, float);
  63 static float jntf(float, float);
  64 static float recurf( float *, float, float * );
  65 float sqrtf(float), sinf(float), cosf(float);
  66 float lgamf(float), expf(float), logf(float), powf(float, float);
  67 float gammaf(float), cbrtf(float), acosf(float);
  68 int airyf(float, float *, float *, float *, float *);
  69 float polevlf(float, float *, int);
  70 #else
  71 float floorf(), j0f(), j1f();
  72 float sqrtf(), sinf(), cosf();
  73 float lgamf(), expf(), logf(), powf(), gammaf();
  74 float cbrtf(), polevlf(), acosf();
  75 void airyf();
  76 static float recurf(), jvsf(), hankelf(), jnxf(), jntf(), jvsf();
  77 #endif
  78
  79
  80 #define fabsf(x) ( (x) < 0 ? -(x) : (x) )
  81
  82 #ifdef ANSIC
  83 float jvf( float nn, float xx )
  84 #else
  85 float jvf( nn, xx )
  86 double nn, xx;
  87 #endif
  88 {
  89 float n, x, k, q, t, y, an, sign;
  90 int i, nint;
  91
  92 n = nn;
  93 x = xx;
  94 nint = 0;       /* Flag for integer n */
  95 sign = 1.0;     /* Flag for sign inversion */
  96 an = fabsf( n );
  97 y = floorf( an );
  98 if( y == an )
  99         {
 100         nint = 1;
 101         i = an - 16384.0 * floorf( an/16384.0 );
 102         if( n < 0.0 )
 103                 {
 104                 if( i & 1 )
 105                         sign = -sign;
 106                 n = an;
 107                 }
 108         if( x < 0.0 )
 109                 {
 110                 if( i & 1 )
 111                         sign = -sign;
 112                 x = -x;
 113                 }
 114         if( n == 0.0 )
 115                 return( j0f(x) );
 116         if( n == 1.0 )
 117                 return( sign * j1f(x) );
 118         }
 119
 120 if( (x < 0.0) && (y != an) )
 121         {
 122         mtherr( "jvf", DOMAIN );
 123         y = 0.0;
 124         goto done;
 125         }
 126
 127 y = fabsf(x);
 128
 129 if( y < MACHEPF )
 130         goto underf;
 131
 132 /* Easy cases - x small compared to n */
 133 t = 3.6 * sqrtf(an);
 134 if( y < t )
 135         return( sign * jvsf(n,x) );
 136
 137 /* x large compared to n */
 138 k = 3.6 * sqrtf(y);
 139 if( (an < k) && (y > 6.0) )
 140         return( sign * hankelf(n,x) );
 141
 142 if( (n > -100) && (n < 14.0) )
 143         {
 144 /* Note: if x is too large, the continued
 145  * fraction will fail; but then the
 146  * Hankel expansion can be used.
 147  */
 148         if( nint != 0 )
 149                 {
 150                 k = 0.0;
 151                 q = recurf( &n, x, &k );
 152                 if( k == 0.0 )
 153                         {
 154                         y = j0f(x)/q;
 155                         goto done;
 156                         }
 157                 if( k == 1.0 )
 158                         {
 159                         y = j1f(x)/q;
 160                         goto done;
 161                         }
 162                 }
 163
 164         if( n >= 0.0 )
 165                 {
 166 /* Recur backwards from a larger value of n
 167  */
 168                 if( y > 1.3 * an )
 169                         goto recurdwn;
 170                 if( an > 1.3 * y )
 171                         goto recurdwn;
 172                 k = n;
 173                 y = 2.0*(y+an+1.0);
 174                 if( (y - n) > 33.0 )
 175                         y = n + 33.0;
 176                 y = n + floorf(y-n);
 177                 q = recurf( &y, x, &k );
 178                 y = jvsf(y,x) * q;
 179                 goto done;
 180                 }
 181 recurdwn:
 182         if( an > (k + 3.0) )
 183                 {
 184 /* Recur backwards from n to k
 185  */
 186                 if( n < 0.0 )
 187                         k = -k;
 188                 q = n - floorf(n);
 189                 k = floorf(k) + q;
 190                 if( n > 0.0 )
 191                         q = recurf( &n, x, &k );
 192                 else
 193                         {
 194                         t = k;
 195                         k = n;
 196                         q = recurf( &t, x, &k );
 197                         k = t;
 198                         }
 199                 if( q == 0.0 )
 200                         {
 201 underf:
 202                         y = 0.0;
 203                         goto done;
 204                         }
 205                 }
 206         else
 207                 {
 208                 k = n;
 209                 q = 1.0;
 210                 }
 211
 212 /* boundary between convergence of
 213  * power series and Hankel expansion
 214  */
 215         t = fabsf(k);
 216         if( t < 26.0 )
 217                 t = (0.0083*t + 0.09)*t + 12.9;
 218         else
 219                 t = 0.9 * t;
 220
 221         if( y > t ) /* y = |x| */
 222                 y = hankelf(k,x);
 223         else
 224                 y = jvsf(k,x);
 225 #if DEBUG
 226 printf( "y = %.16e, q = %.16e\n", y, q );
 227 #endif
 228         if( n > 0.0 )
 229                 y /= q;
 230         else
 231                 y *= q;
 232         }
 233
 234 else
 235         {
 236 /* For large positive n, use the uniform expansion
 237  * or the transitional expansion.
 238  * But if x is of the order of n**2,
 239  * these may blow up, whereas the
 240  * Hankel expansion will then work.
 241  */
 242         if( n < 0.0 )
 243                 {
 244                 mtherr( "jvf", TLOSS );
 245                 y = 0.0;
 246                 goto done;
 247                 }
 248         t = y/an;
 249         t /= an;
 250         if( t > 0.3 )
 251                 y = hankelf(n,x);
 252         else
 253                 y = jnxf(n,x);
 254         }
 255
 256 done:   return( sign * y);
 257 }
 258 \f
 259 /* Reduce the order by backward recurrence.
 260  * AMS55 #9.1.27 and 9.1.73.
 261  */
 262
 263 #ifdef ANSIC
 264 static float recurf( float *n, float xx, float *newn )
 265 #else
 266 static float recurf( n, xx, newn )
 267 float *n;
 268 double xx;
 269 float *newn;
 270 #endif
 271 {
 272 float x, pkm2, pkm1, pk, pkp1, qkm2, qkm1;
 273 float k, ans, qk, xk, yk, r, t, kf, xinv;
 274 static float big = BIG;
 275 int nflag, ctr;
 276
 277 x = xx;
 278 /* continued fraction for Jn(x)/Jn-1(x)  */
 279 if( *n < 0.0 )
 280         nflag = 1;
 281 else
 282         nflag = 0;
 283
 284 fstart:
 285
 286 #if DEBUG
 287 printf( "n = %.6e, newn = %.6e, cfrac = ", *n, *newn );
 288 #endif
 289
 290 pkm2 = 0.0;
 291 qkm2 = 1.0;
 292 pkm1 = x;
 293 qkm1 = *n + *n;
 294 xk = -x * x;
 295 yk = qkm1;
 296 ans = 1.0;
 297 ctr = 0;
 298 do
 299         {
 300         yk += 2.0;
 301         pk = pkm1 * yk +  pkm2 * xk;
 302         qk = qkm1 * yk +  qkm2 * xk;
 303         pkm2 = pkm1;
 304         pkm1 = pk;
 305         qkm2 = qkm1;
 306         qkm1 = qk;
 307         if( qk != 0 )
 308                 r = pk/qk;
 309         else
 310                 r = 0.0;
 311         if( r != 0 )
 312                 {
 313                 t = fabsf( (ans - r)/r );
 314                 ans = r;
 315                 }
 316         else
 317                 t = 1.0;
 318
 319         if( t < MACHEPF )
 320                 goto done;
 321
 322         if( fabsf(pk) > big )
 323                 {
 324                 pkm2 *= MACHEPF;
 325                 pkm1 *= MACHEPF;
 326                 qkm2 *= MACHEPF;
 327                 qkm1 *= MACHEPF;
 328                 }
 329         }
 330 while( t > MACHEPF );
 331
 332 done:
 333
 334 #if DEBUG
 335 printf( "%.6e\n", ans );
 336 #endif
 337
 338 /* Change n to n-1 if n < 0 and the continued fraction is small
 339  */
 340 if( nflag > 0 )
 341         {
 342         if( fabsf(ans) < 0.125 )
 343                 {
 344                 nflag = -1;
 345                 *n = *n - 1.0;
 346                 goto fstart;
 347                 }
 348         }
 349
 350
 351 kf = *newn;
 352
 353 /* backward recurrence
 354  *              2k
 355  *  J   (x)  =  --- J (x)  -  J   (x)
 356  *   k-1         x   k         k+1
 357  */
 358
 359 pk = 1.0;
 360 pkm1 = 1.0/ans;
 361 k = *n - 1.0;
 362 r = 2 * k;
 363 xinv = 1.0/x;
 364 do
 365         {
 366         pkm2 = (pkm1 * r  -  pk * x) * xinv;
 367         pkp1 = pk;
 368         pk = pkm1;
 369         pkm1 = pkm2;
 370         r -= 2.0;
 371 #if 0
 372         t = fabsf(pkp1) + fabsf(pk);
 373         if( (k > (kf + 2.5)) && (fabsf(pkm1) < 0.25*t) )
 374                 {
 375                 k -= 1.0;
 376                 t = x*x;
 377                 pkm2 = ( (r*(r+2.0)-t)*pk - r*x*pkp1 )/t;
 378                 pkp1 = pk;
 379                 pk = pkm1;
 380                 pkm1 = pkm2;
 381                 r -= 2.0;
 382                 }
 383 #endif
 384         k -= 1.0;
 385         }
 386 while( k > (kf + 0.5) );
 387
 388 #if 0
 389 /* Take the larger of the last two iterates
 390  * on the theory that it may have less cancellation error.
 391  */
 392 if( (kf >= 0.0) && (fabsf(pk) > fabsf(pkm1)) )
 393         {
 394         k += 1.0;
 395         pkm2 = pk;
 396         }
 397 #endif
 398
 399 *newn = k;
 400 #if DEBUG
 401 printf( "newn %.6e\n", k );
 402 #endif
 403 return( pkm2 );
 404 }
 405
 406
 407
 408 /* Ascending power series for Jv(x).
 409  * AMS55 #9.1.10.
 410  */
 411
 412 #ifdef ANSIC
 413 static float jvsf( float nn, float xx )
 414 #else
 415 static float jvsf( nn, xx )
 416 double nn, xx;
 417 #endif
 418 {
 419 float n, x, t, u, y, z, k, ay;
 420
 421 #if DEBUG
 422 printf( "jvsf: " );
 423 #endif
 424 n = nn;
 425 x = xx;
 426 z = -0.25 * x * x;
 427 u = 1.0;
 428 y = u;
 429 k = 1.0;
 430 t = 1.0;
 431
 432 while( t > MACHEPF )
 433         {
 434         u *= z / (k * (n+k));
 435         y += u;
 436         k += 1.0;
 437         t = fabsf(u);
 438         if( (ay = fabsf(y)) > 1.0 )
 439                 t /= ay;
 440         }
 441
 442 if( x < 0.0 )
 443         {
 444         y = y * powf( 0.5 * x, n ) / gammaf( n + 1.0 );
 445         }
 446 else
 447         {
 448         t = n * logf(0.5*x) - lgamf(n + 1.0);
 449         if( t < -MAXLOGF )
 450                 {
 451                 return( 0.0 );
 452                 }
 453         if( t > MAXLOGF )
 454                 {
 455                 t = logf(y) + t;
 456                 if( t > MAXLOGF )
 457                         {
 458                         mtherr( "jvf", OVERFLOW );
 459                         return( MAXNUMF );
 460                         }
 461                 else
 462                         {
 463                         y = sgngamf * expf(t);
 464                         return(y);
 465                         }
 466                 }
 467         y = sgngamf * y * expf( t );
 468         }
 469 #if DEBUG
 470 printf( "y = %.8e\n", y );
 471 #endif
 472 return(y);
 473 }
 474 \f
 475 /* Hankel's asymptotic expansion
 476  * for large x.
 477  * AMS55 #9.2.5.
 478  */
 479 #ifdef ANSIC
 480 static float hankelf( float nn, float xx )
 481 #else
 482 static float hankelf( nn, xx )
 483 double nn, xx;
 484 #endif
 485 {
 486 float n, x, t, u, z, k, sign, conv;
 487 float p, q, j, m, pp, qq;
 488 int flag;
 489
 490 #if DEBUG
 491 printf( "hankelf: " );
 492 #endif
 493 n = nn;
 494 x = xx;
 495 m = 4.0*n*n;
 496 j = 1.0;
 497 z = 8.0 * x;
 498 k = 1.0;
 499 p = 1.0;
 500 u = (m - 1.0)/z;
 501 q = u;
 502 sign = 1.0;
 503 conv = 1.0;
 504 flag = 0;
 505 t = 1.0;
 506 pp = 1.0e38;
 507 qq = 1.0e38;
 508
 509 while( t > MACHEPF )
 510         {
 511         k += 2.0;
 512         j += 1.0;
 513         sign = -sign;
 514         u *= (m - k * k)/(j * z);
 515         p += sign * u;
 516         k += 2.0;
 517         j += 1.0;
 518         u *= (m - k * k)/(j * z);
 519         q += sign * u;
 520         t = fabsf(u/p);
 521         if( t < conv )
 522                 {
 523                 conv = t;
 524                 qq = q;
 525                 pp = p;
 526                 flag = 1;
 527                 }
 528 /* stop if the terms start getting larger */
 529         if( (flag != 0) && (t > conv) )
 530                 {
 531 #if DEBUG
 532                 printf( "Hankel: convergence to %.4E\n", conv );
 533 #endif
 534                 goto hank1;
 535                 }
 536         }
 537
 538 hank1:
 539 u = x - (0.5*n + 0.25) * PIF;
 540 t = sqrtf( 2.0/(PIF*x) ) * ( pp * cosf(u) - qq * sinf(u) );
 541 return( t );
 542 }
 543 \f
 544
 545 /* Asymptotic expansion for large n.
 546  * AMS55 #9.3.35.
 547  */
 548
 549 static float lambda[] = {
 550   1.0,
 551   1.041666666666666666666667E-1,
 552   8.355034722222222222222222E-2,
 553   1.282265745563271604938272E-1,
 554   2.918490264641404642489712E-1,
 555   8.816272674437576524187671E-1,
 556   3.321408281862767544702647E+0,
 557   1.499576298686255465867237E+1,
 558   7.892301301158651813848139E+1,
 559   4.744515388682643231611949E+2,
 560   3.207490090890661934704328E+3
 561 };
 562 static float mu[] = {
 563   1.0,
 564  -1.458333333333333333333333E-1,
 565  -9.874131944444444444444444E-2,
 566  -1.433120539158950617283951E-1,
 567  -3.172272026784135480967078E-1,
 568  -9.424291479571202491373028E-1,
 569  -3.511203040826354261542798E+0,
 570  -1.572726362036804512982712E+1,
 571  -8.228143909718594444224656E+1,
 572  -4.923553705236705240352022E+2,
 573  -3.316218568547972508762102E+3
 574 };
 575 static float P1[] = {
 576  -2.083333333333333333333333E-1,
 577   1.250000000000000000000000E-1
 578 };
 579 static float P2[] = {
 580   3.342013888888888888888889E-1,
 581  -4.010416666666666666666667E-1,
 582   7.031250000000000000000000E-2
 583 };
 584 static float P3[] = {
 585  -1.025812596450617283950617E+0,
 586   1.846462673611111111111111E+0,
 587  -8.912109375000000000000000E-1,
 588   7.324218750000000000000000E-2
 589 };
 590 static float P4[] = {
 591   4.669584423426247427983539E+0,
 592  -1.120700261622299382716049E+1,
 593   8.789123535156250000000000E+0,
 594  -2.364086914062500000000000E+0,
 595   1.121520996093750000000000E-1
 596 };
 597 static float P5[] = {
 598  -2.8212072558200244877E1,
 599   8.4636217674600734632E1,
 600  -9.1818241543240017361E1,
 601   4.2534998745388454861E1,
 602  -7.3687943594796316964E0,
 603   2.27108001708984375E-1
 604 };
 605 static float P6[] = {
 606   2.1257013003921712286E2,
 607  -7.6525246814118164230E2,
 608   1.0599904525279998779E3,
 609  -6.9957962737613254123E2,
 610   2.1819051174421159048E2,
 611  -2.6491430486951555525E1,
 612   5.7250142097473144531E-1
 613 };
 614 static float P7[] = {
 615  -1.9194576623184069963E3,
 616   8.0617221817373093845E3,
 617  -1.3586550006434137439E4,
 618   1.1655393336864533248E4,
 619  -5.3056469786134031084E3,
 620   1.2009029132163524628E3,
 621  -1.0809091978839465550E2,
 622   1.7277275025844573975E0
 623 };
 624
 625
 626 #ifdef ANSIC
 627 static float jnxf( float nn, float xx )
 628 #else
 629 static float jnxf( nn, xx )
 630 double nn, xx;
 631 #endif
 632 {
 633 float n, x, zeta, sqz, zz, zp, np;
 634 float cbn, n23, t, z, sz;
 635 float pp, qq, z32i, zzi;
 636 float ak, bk, akl, bkl;
 637 int sign, doa, dob, nflg, k, s, tk, tkp1, m;
 638 static float u[8];
 639 static float ai, aip, bi, bip;
 640
 641 n = nn;
 642 x = xx;
 643 /* Test for x very close to n.
 644  * Use expansion for transition region if so.
 645  */
 646 cbn = cbrtf(n);
 647 z = (x - n)/cbn;
 648 if( (fabsf(z) <= 0.7) || (n < 0.0) )
 649         return( jntf(n,x) );
 650 z = x/n;
 651 zz = 1.0 - z*z;
 652 if( zz == 0.0 )
 653         return(0.0);
 654
 655 if( zz > 0.0 )
 656         {
 657         sz = sqrtf( zz );
 658         t = 1.5 * (logf( (1.0+sz)/z ) - sz );   /* zeta ** 3/2          */
 659         zeta = cbrtf( t * t );
 660         nflg = 1;
 661         }
 662 else
 663         {
 664         sz = sqrtf(-zz);
 665         t = 1.5 * (sz - acosf(1.0/z));
 666         zeta = -cbrtf( t * t );
 667         nflg = -1;
 668         }
 669 z32i = fabsf(1.0/t);
 670 sqz = cbrtf(t);
 671
 672 /* Airy function */
 673 n23 = cbrtf( n * n );
 674 t = n23 * zeta;
 675
 676 #if DEBUG
 677 printf("zeta %.5E, Airyf(%.5E)\n", zeta, t );
 678 #endif
 679 airyf( t, &ai, &aip, &bi, &bip );
 680
 681 /* polynomials in expansion */
 682 u[0] = 1.0;
 683 zzi = 1.0/zz;
 684 u[1] = polevlf( zzi, P1, 1 )/sz;
 685 u[2] = polevlf( zzi, P2, 2 )/zz;
 686 u[3] = polevlf( zzi, P3, 3 )/(sz*zz);
 687 pp = zz*zz;
 688 u[4] = polevlf( zzi, P4, 4 )/pp;
 689 u[5] = polevlf( zzi, P5, 5 )/(pp*sz);
 690 pp *= zz;
 691 u[6] = polevlf( zzi, P6, 6 )/pp;
 692 u[7] = polevlf( zzi, P7, 7 )/(pp*sz);
 693
 694 #if DEBUG
 695 for( k=0; k<=7; k++ )
 696         printf( "u[%d] = %.5E\n", k, u[k] );
 697 #endif
 698
 699 pp = 0.0;
 700 qq = 0.0;
 701 np = 1.0;
 702 /* flags to stop when terms get larger */
 703 doa = 1;
 704 dob = 1;
 705 akl = MAXNUMF;
 706 bkl = MAXNUMF;
 707
 708 for( k=0; k<=3; k++ )
 709         {
 710         tk = 2 * k;
 711         tkp1 = tk + 1;
 712         zp = 1.0;
 713         ak = 0.0;
 714         bk = 0.0;
 715         for( s=0; s<=tk; s++ )
 716                 {
 717                 if( doa )
 718                         {
 719                         if( (s & 3) > 1 )
 720                                 sign = nflg;
 721                         else
 722                                 sign = 1;
 723                         ak += sign * mu[s] * zp * u[tk-s];
 724                         }
 725
 726                 if( dob )
 727                         {
 728                         m = tkp1 - s;
 729                         if( ((m+1) & 3) > 1 )
 730                                 sign = nflg;
 731                         else
 732                                 sign = 1;
 733                         bk += sign * lambda[s] * zp * u[m];
 734                         }
 735                 zp *= z32i;
 736                 }
 737
 738         if( doa )
 739                 {
 740                 ak *= np;
 741                 t = fabsf(ak);
 742                 if( t < akl )
 743                         {
 744                         akl = t;
 745                         pp += ak;
 746                         }
 747                 else
 748                         doa = 0;
 749                 }
 750
 751         if( dob )
 752                 {
 753                 bk += lambda[tkp1] * zp * u[0];
 754                 bk *= -np/sqz;
 755                 t = fabsf(bk);
 756                 if( t < bkl )
 757                         {
 758                         bkl = t;
 759                         qq += bk;
 760                         }
 761                 else
 762                         dob = 0;
 763                 }
 764 #if DEBUG
 765         printf("a[%d] %.5E, b[%d] %.5E\n", k, ak, k, bk );
 766 #endif
 767         if( np < MACHEPF )
 768                 break;
 769         np /= n*n;
 770         }
 771
 772 /* normalizing factor ( 4*zeta/(1 - z**2) )**1/4        */
 773 t = 4.0 * zeta/zz;
 774 t = sqrtf( sqrtf(t) );
 775
 776 t *= ai*pp/cbrtf(n)  +  aip*qq/(n23*n);
 777 return(t);
 778 }
 779 \f
 780 /* Asymptotic expansion for transition region,
 781  * n large and x close to n.
 782  * AMS55 #9.3.23.
 783  */
 784
 785 static float PF2[] = {
 786  -9.0000000000000000000e-2,
 787   8.5714285714285714286e-2
 788 };
 789 static float PF3[] = {
 790   1.3671428571428571429e-1,
 791  -5.4920634920634920635e-2,
 792  -4.4444444444444444444e-3
 793 };
 794 static float PF4[] = {
 795   1.3500000000000000000e-3,
 796  -1.6036054421768707483e-1,
 797   4.2590187590187590188e-2,
 798   2.7330447330447330447e-3
 799 };
 800 static float PG1[] = {
 801  -2.4285714285714285714e-1,
 802   1.4285714285714285714e-2
 803 };
 804 static float PG2[] = {
 805  -9.0000000000000000000e-3,
 806   1.9396825396825396825e-1,
 807  -1.1746031746031746032e-2
 808 };
 809 static float PG3[] = {
 810   1.9607142857142857143e-2,
 811  -1.5983694083694083694e-1,
 812   6.3838383838383838384e-3
 813 };
 814
 815
 816 #ifdef ANSIC
 817 static float jntf( float nn, float xx )
 818 #else
 819 static float jntf( nn, xx )
 820 double nn, xx;
 821 #endif
 822 {
 823 float n, x, z, zz, z3;
 824 float cbn, n23, cbtwo;
 825 float ai, aip, bi, bip; /* Airy functions */
 826 float nk, fk, gk, pp, qq;
 827 float F[5], G[4];
 828 int k;
 829
 830 n = nn;
 831 x = xx;
 832 cbn = cbrtf(n);
 833 z = (x - n)/cbn;
 834 cbtwo = cbrtf( 2.0 );
 835
 836 /* Airy function */
 837 zz = -cbtwo * z;
 838 airyf( zz, &ai, &aip, &bi, &bip );
 839
 840 /* polynomials in expansion */
 841 zz = z * z;
 842 z3 = zz * z;
 843 F[0] = 1.0;
 844 F[1] = -z/5.0;
 845 F[2] = polevlf( z3, PF2, 1 ) * zz;
 846 F[3] = polevlf( z3, PF3, 2 );
 847 F[4] = polevlf( z3, PF4, 3 ) * z;
 848 G[0] = 0.3 * zz;
 849 G[1] = polevlf( z3, PG1, 1 );
 850 G[2] = polevlf( z3, PG2, 2 ) * z;
 851 G[3] = polevlf( z3, PG3, 2 ) * zz;
 852 #if DEBUG
 853 for( k=0; k<=4; k++ )
 854         printf( "F[%d] = %.5E\n", k, F[k] );
 855 for( k=0; k<=3; k++ )
 856         printf( "G[%d] = %.5E\n", k, G[k] );
 857 #endif
 858 pp = 0.0;
 859 qq = 0.0;
 860 nk = 1.0;
 861 n23 = cbrtf( n * n );
 862
 863 for( k=0; k<=4; k++ )
 864         {
 865         fk = F[k]*nk;
 866         pp += fk;
 867         if( k != 4 )
 868                 {
 869                 gk = G[k]*nk;
 870                 qq += gk;
 871                 }
 872 #if DEBUG
 873         printf("fk[%d] %.5E, gk[%d] %.5E\n", k, fk, k, gk );
 874 #endif
 875         nk /= n23;
 876         }
 877
 878 fk = cbtwo * ai * pp/cbn  +  cbrtf(4.0) * aip * qq/n;
 879 return(fk);
 880 }