Return to jn.c CVS log

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

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

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

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