Return to jn.c CVS log

Up to [CSRG BSD Unix] / 43BSD / usr.lib / libm

Annotation of 43BSD/usr.lib/libm/jn.c, revision 1.1.1.1

1.1       root        1: #ifndef lint
                      2: static char sccsid[] = "@(#)jn.c       4.2 (Berkeley) 8/21/85";
                      3: #endif not lint
                      4: 
                      5: /*
                      6:        floating point Bessel's function of
                      7:        the first and second kinds and of
                      8:        integer order.
                      9: 
                     10:        int n;
                     11:        double x;
                     12:        jn(n,x);
                     13: 
                     14:        returns the value of Jn(x) for all
                     15:        integer values of n and all real values
                     16:        of x.
                     17: 
                     18:        There are no error returns.
                     19:        Calls j0, j1.
                     20: 
                     21:        For n=0, j0(x) is called,
                     22:        for n=1, j1(x) is called,
                     23:        for n<x, forward recursion us used starting
                     24:        from values of j0(x) and j1(x).
                     25:        for n>x, a continued fraction approximation to
                     26:        j(n,x)/j(n-1,x) is evaluated and then backward
                     27:        recursion is used starting from a supposed value
                     28:        for j(n,x). The resulting value of j(0,x) is
                     29:        compared with the actual value to correct the
                     30:        supposed value of j(n,x).
                     31: 
                     32:        yn(n,x) is similar in all respects, except
                     33:        that forward recursion is used for all
                     34:        values of n>1.
                     35: */
                     36: 
                     37: #include <math.h>
                     38: #ifdef VAX
                     39: #include <errno.h>
                     40: #else  /* IEEE double */
                     41: static double zero = 0.e0;
                     42: #endif
                     43: 
                     44: double
                     45: jn(n,x) int n; double x;{
                     46:        int i;
                     47:        double a, b, temp;
                     48:        double xsq, t;
                     49:        double j0(), j1();
                     50: 
                     51:        if(n<0){
                     52:                n = -n;
                     53:                x = -x;
                     54:        }
                     55:        if(n==0) return(j0(x));
                     56:        if(n==1) return(j1(x));
                     57:        if(x == 0.) return(0.);
                     58:        if(n>x) goto recurs;
                     59: 
                     60:        a = j0(x);
                     61:        b = j1(x);
                     62:        for(i=1;i<n;i++){
                     63:                temp = b;
                     64:                b = (2.*i/x)*b - a;
                     65:                a = temp;
                     66:        }
                     67:        return(b);
                     68: 
                     69: recurs:
                     70:        xsq = x*x;
                     71:        for(t=0,i=n+16;i>n;i--){
                     72:                t = xsq/(2.*i - t);
                     73:        }
                     74:        t = x/(2.*n-t);
                     75: 
                     76:        a = t;
                     77:        b = 1;
                     78:        for(i=n-1;i>0;i--){
                     79:                temp = b;
                     80:                b = (2.*i/x)*b - a;
                     81:                a = temp;
                     82:        }
                     83:        return(t*j0(x)/b);
                     84: }
                     85: 
                     86: double
                     87: yn(n,x) int n; double x;{
                     88:        int i;
                     89:        int sign;
                     90:        double a, b, temp;
                     91:        double y0(), y1();
                     92: 
                     93:        if (x <= 0) {
                     94: #ifdef VAX
                     95:                extern double infnan();
                     96:                return(infnan(EDOM));   /* NaN */
                     97: #else  /* IEEE double */
                     98:                return(zero/zero);      /* IEEE machines: invalid operation */
                     99: #endif
                    100:        }
                    101:        sign = 1;
                    102:        if(n<0){
                    103:                n = -n;
                    104:                if(n%2 == 1) sign = -1;
                    105:        }
                    106:        if(n==0) return(y0(x));
                    107:        if(n==1) return(sign*y1(x));
                    108: 
                    109:        a = y0(x);
                    110:        b = y1(x);
                    111:        for(i=1;i<n;i++){
                    112:                temp = b;
                    113:                b = (2.*i/x)*b - a;
                    114:                a = temp;
                    115:        }
                    116:        return(sign*b);
                    117: }

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