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