Return to jn.c CVS log

Up to [CSRG BSD Unix] / 43BSDReno / lib / libm / common_source

Annotation of 43BSDReno/lib/libm/common_source/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.3 (Berkeley) 9/22/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 "mathimpl.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: 
                     56:        if(n<0){
                     57:                n = -n;
                     58:                x = -x;
                     59:        }
                     60:        if(n==0) return(j0(x));
                     61:        if(n==1) return(j1(x));
                     62:        if(x == 0.) return(0.);
                     63:        if(n>x) goto recurs;
                     64: 
                     65:        a = j0(x);
                     66:        b = j1(x);
                     67:        for(i=1;i<n;i++){
                     68:                temp = b;
                     69:                b = (2.*i/x)*b - a;
                     70:                a = temp;
                     71:        }
                     72:        return(b);
                     73: 
                     74: recurs:
                     75:        xsq = x*x;
                     76:        for(t=0,i=n+16;i>n;i--){
                     77:                t = xsq/(2.*i - t);
                     78:        }
                     79:        t = x/(2.*n-t);
                     80: 
                     81:        a = t;
                     82:        b = 1;
                     83:        for(i=n-1;i>0;i--){
                     84:                temp = b;
                     85:                b = (2.*i/x)*b - a;
                     86:                a = temp;
                     87:        }
                     88:        return(t*j0(x)/b);
                     89: }
                     90: 
                     91: double
                     92: yn(n,x) int n; double x;{
                     93:        int i;
                     94:        int sign;
                     95:        double a, b, temp;
                     96: 
                     97:        if (x <= 0) {
                     98: #if defined(vax)||defined(tahoe)
                     99:                return(infnan(EDOM));   /* NaN */
                    100: #else  /* defined(vax)||defined(tahoe) */
                    101:                return(zero/zero);      /* IEEE machines: invalid operation */
                    102: #endif /* defined(vax)||defined(tahoe) */
                    103:        }
                    104:        sign = 1;
                    105:        if(n<0){
                    106:                n = -n;
                    107:                if(n%2 == 1) sign = -1;
                    108:        }
                    109:        if(n==0) return(y0(x));
                    110:        if(n==1) return(sign*y1(x));
                    111: 
                    112:        a = y0(x);
                    113:        b = y1(x);
                    114:        for(i=1;i<n;i++){
                    115:                temp = b;
                    116:                b = (2.*i/x)*b - a;
                    117:                a = temp;
                    118:        }
                    119:        return(sign*b);
                    120: }

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