home *** CD-ROM | disk | FTP | other *** search
/ Il CD di internet / CD.iso / SOURCE / D / LIBC / LIBC-4.6 / LIBC-4 / libc-linux / sysdeps / linux / i386 / math / jn.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-19  |  1.9 KB  |  120 lines
  1. #include <math.h>
  2.  
  3. /*
  4.     floating point Bessel's function of
  5.     the first and second kinds and of
  6.     integer order.
  7.  
  8.     int n;
  9.     double x;
  10.     jn(n,x);
  11.  
  12.     returns the value of Jn(x) for all
  13.     integer values of n and all real values
  14.     of x.
  15.  
  16.     There are no error returns.
  17.     Calls j0, j1.
  18.  
  19.     For n=0, j0(x) is called,
  20.     for n=1, j1(x) is called,
  21.     for n<x, forward recursion us used starting
  22.     from values of j0(x) and j1(x).
  23.     for n>x, a continued fraction approximation to
  24.     j(n,x)/j(n-1,x) is evaluated and then backward
  25.     recursion is used starting from a supposed value
  26.     for j(n,x). The resulting value of j(0,x) is
  27.     compared with the actual value to correct the
  28.     supposed value of j(n,x).
  29.  
  30.     yn(n,x) is similar in all respects, except
  31.     that forward recursion is used for all
  32.     values of n>1.
  33. */
  34.  
  35. double jn(int n, double x)
  36. {
  37.     int i, sign;
  38.     double a, b, temp;
  39.     double xsq, t;
  40.  
  41. /*
  42.                   n
  43.    J  (x)  =  (-1)   J (x)
  44.     -n                n
  45.                   n
  46.    J (-x)  =  (-1)   J (x)
  47.     n                 n
  48.  */
  49.     sign = 1;
  50.     if(n<0){
  51.         n = -n;
  52.         if( n & 1 )
  53.           sign = -1;
  54.     }
  55.     if( x < 0.0 )
  56.       {
  57.         if( n & 1 )
  58.           sign = -sign;
  59.         x = -x;
  60.       }
  61.     if(n==0) return(j0(x));
  62.     if(n==1) return(sign*j1(x));
  63.     if(x == 0.) return(0.);
  64.     if(n>x) goto recurs;
  65.  
  66. /* Note, forward recurrence for Jn is numerically unstable.
  67.    This should be fixed sometime.  */
  68.     a = j0(x);
  69.     b = j1(x);
  70.     for(i=1;i<n;i++){
  71.         temp = b;
  72.         b = (2.*i/x)*b - a;
  73.         a = temp;
  74.     }
  75.     return(sign*b);
  76.  
  77. recurs:
  78.     xsq = x*x;
  79.     for(t=0,i=n+16;i>n;i--){
  80.         t = xsq/(2.*i - t);
  81.     }
  82.     t = x/(2.*n-t);
  83.  
  84.     a = t;
  85.     b = 1;
  86.     for(i=n-1;i>0;i--){
  87.         temp = b;
  88.         b = (2.*i/x)*b - a;
  89.         a = temp;
  90.     }
  91.     return(sign*t*j0(x)/b);
  92. }
  93.  
  94. double yn(int n, double x)
  95. {
  96.     int i;
  97.     int sign;
  98.     double a, b, temp;
  99.  
  100.     if (x <= 0) {
  101.         return(-HUGE_VAL);
  102.     }
  103.     sign = 1;
  104.     if(n<0){
  105.         n = -n;
  106.         if(n%2 == 1) sign = -1;
  107.     }
  108.     if(n==0) return(y0(x));
  109.     if(n==1) return(sign*y1(x));
  110.  
  111.     a = y0(x);
  112.     b = y1(x);
  113.     for(i=1;i<n;i++){
  114.         temp = b;
  115.         b = (2.*i/x)*b - a;
  116.         a = temp;
  117.     }
  118.     return(sign*b);
  119. }
  120.