Return to besjn.c CVS log

Up to [Research Unix] / researchv9 / libc / math

Annotation of researchv9/libc/math/besjn.c, revision 1.1.1.1

1.1       root        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 j0, j1.
                     16: 
                     17:        For n=0, j0(x) is called,
                     18:        for n=1, j1(x) is called,
                     19:        for n<x, forward recursion us used starting
                     20:        from values of j0(x) and j1(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 <math.h>
                     34: #include <errno.h>
                     35: 
                     36: int    errno;
                     37: 
                     38: double
                     39: besjn(n,x) int n; double x;{
                     40:        int i;
                     41:        double a, b, temp;
                     42:        double xsq, t;
                     43:        double j0(), j1();
                     44: 
                     45:        if(n<0){
                     46:                n = -n;
                     47:                x = -x;
                     48:        }
                     49:        if(n==0) return(j0(x));
                     50:        if(n==1) return(j1(x));
                     51:        if(x == 0.) return(0.);
                     52:        if(n>x) goto recurs;
                     53: 
                     54:        a = j0(x);
                     55:        b = j1(x);
                     56:        for(i=1;i<n;i++){
                     57:                temp = b;
                     58:                b = (2.*i/x)*b - a;
                     59:                a = temp;
                     60:        }
                     61:        return(b);
                     62: 
                     63: recurs:
                     64:        xsq = x*x;
                     65:        for(t=0,i=n+16;i>n;i--){
                     66:                t = xsq/(2.*i - t);
                     67:        }
                     68:        t = x/(2.*n-t);
                     69: 
                     70:        a = t;
                     71:        b = 1;
                     72:        for(i=n-1;i>0;i--){
                     73:                temp = b;
                     74:                b = (2.*i/x)*b - a;
                     75:                a = temp;
                     76:        }
                     77:        return(t*j0(x)/b);
                     78: }
                     79: 
                     80: double
                     81: besyn(n,x) int n; double x;{
                     82:        int i;
                     83:        int sign;
                     84:        double a, b, temp;
                     85:        double y0(), y1();
                     86: 
                     87:        if (x <= 0) {
                     88:                errno = EDOM;
                     89:                return(-HUGE);
                     90:        }
                     91:        sign = 1;
                     92:        if(n<0){
                     93:                n = -n;
                     94:                if(n%2 == 1) sign = -1;
                     95:        }
                     96:        if(n==0) return(y0(x));
                     97:        if(n==1) return(sign*y1(x));
                     98: 
                     99:        a = y0(x);
                    100:        b = y1(x);
                    101:        for(i=1;i<n;i++){
                    102:                temp = b;
                    103:                b = (2.*i/x)*b - a;
                    104:                a = temp;
                    105:        }
                    106:        return(sign*b);
                    107: }

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.