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     ! 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.