libm/ldouble/jnl.c

   1 /*                                                      jnl.c
   2  *
   3  *      Bessel function of integer order
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * int n;
  10  * long double x, y, jnl();
  11  *
  12  * y = jnl( n, x );
  13  *
  14  *
  15  *
  16  * DESCRIPTION:
  17  *
  18  * Returns Bessel function of order n, where n is a
  19  * (possibly negative) integer.
  20  *
  21  * The ratio of jn(x) to j0(x) is computed by backward
  22  * recurrence.  First the ratio jn/jn-1 is found by a
  23  * continued fraction expansion.  Then the recurrence
  24  * relating successive orders is applied until j0 or j1 is
  25  * reached.
  26  *
  27  * If n = 0 or 1 the routine for j0 or j1 is called
  28  * directly.
  29  *
  30  *
  31  *
  32  * ACCURACY:
  33  *
  34  *                      Absolute error:
  35  * arithmetic   domain      # trials      peak         rms
  36  *    IEEE     -30, 30        5000       3.3e-19     4.7e-20
  37  *
  38  *
  39  * Not suitable for large n or x.
  40  *
  41  */
  42 \f
  43 /*                                                      jn.c
  44 Cephes Math Library Release 2.0:  April, 1987
  45 Copyright 1984, 1987 by Stephen L. Moshier
  46 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  47 */
  48 #include <math.h>
  49
  50 extern long double MACHEPL;
  51 #ifdef ANSIPROT
  52 extern long double fabsl ( long double );
  53 extern long double j0l ( long double );
  54 extern long double j1l ( long double );
  55 #else
  56 long double fabsl(), j0l(), j1l();
  57 #endif
  58
  59 long double jnl( n, x )
  60 int n;
  61 long double x;
  62 {
  63 long double pkm2, pkm1, pk, xk, r, ans;
  64 int k, sign;
  65
  66 if( n < 0 )
  67         {
  68         n = -n;
  69         if( (n & 1) == 0 )      /* -1**n */
  70                 sign = 1;
  71         else
  72                 sign = -1;
  73         }
  74 else
  75         sign = 1;
  76
  77 if( x < 0.0L )
  78         {
  79         if( n & 1 )
  80                 sign = -sign;
  81         x = -x;
  82         }
  83
  84
  85 if( n == 0 )
  86         return( sign * j0l(x) );
  87 if( n == 1 )
  88         return( sign * j1l(x) );
  89 if( n == 2 )
  90         return( sign * (2.0L * j1l(x) / x  -  j0l(x)) );
  91
  92 if( x < MACHEPL )
  93         return( 0.0L );
  94
  95 /* continued fraction */
  96 k = 53;
  97 pk = 2 * (n + k);
  98 ans = pk;
  99 xk = x * x;
 100
 101 do
 102         {
 103         pk -= 2.0L;
 104         ans = pk - (xk/ans);
 105         }
 106 while( --k > 0 );
 107 ans = x/ans;
 108
 109 /* backward recurrence */
 110
 111 pk = 1.0L;
 112 pkm1 = 1.0L/ans;
 113 k = n-1;
 114 r = 2 * k;
 115
 116 do
 117         {
 118         pkm2 = (pkm1 * r  -  pk * x) / x;
 119         pk = pkm1;
 120         pkm1 = pkm2;
 121         r -= 2.0L;
 122         }
 123 while( --k > 0 );
 124
 125 if( fabsl(pk) > fabsl(pkm1) )
 126         ans = j1l(x)/pk;
 127 else
 128         ans = j0l(x)/pkm1;
 129 return( sign * ans );
 130 }