Return to jn.c CVS log

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

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

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

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