Return to j1.c CVS log

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

Annotation of 43BSDTahoe/usr.lib/libm/j1.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[] = "@(#)j1.c       5.2 (Berkeley) 4/29/88";
        !             9: #endif /* not lint */
        !            10: 
        !            11: /*
        !            12:        floating point Bessel's function
        !            13:        of the first and second kinds
        !            14:        of order one
        !            15: 
        !            16:        j1(x) returns the value of J1(x)
        !            17:        for all real values of x.
        !            18: 
        !            19:        There are no error returns.
        !            20:        Calls sin, cos, sqrt.
        !            21: 
        !            22:        There is a niggling bug in J1 which
        !            23:        causes errors up to 2e-16 for x in the
        !            24:        interval [-8,8].
        !            25:        The bug is caused by an inappropriate order
        !            26:        of summation of the series.  rhm will fix it
        !            27:        someday.
        !            28: 
        !            29:        Coefficients are from Hart & Cheney.
        !            30:        #6050 (20.98D)
        !            31:        #6750 (19.19D)
        !            32:        #7150 (19.35D)
        !            33: 
        !            34:        y1(x) returns the value of Y1(x)
        !            35:        for positive real values of x.
        !            36:        For x<=0, if on the VAX, error number EDOM is set and
        !            37:        the reserved operand fault is generated;
        !            38:        otherwise (an IEEE machine) an invalid operation is performed.
        !            39: 
        !            40:        Calls sin, cos, sqrt, log, j1.
        !            41: 
        !            42:        The values of Y1 have not been checked
        !            43:        to more than ten places.
        !            44: 
        !            45:        Coefficients are from Hart & Cheney.
        !            46:        #6447 (22.18D)
        !            47:        #6750 (19.19D)
        !            48:        #7150 (19.35D)
        !            49: */
        !            50: 
        !            51: #include <math.h>
        !            52: #if defined(vax)||defined(tahoe)
        !            53: #include <errno.h>
        !            54: #else  /* defined(vax)||defined(tahoe) */
        !            55: static double zero = 0.e0;
        !            56: #endif /* defined(vax)||defined(tahoe) */
        !            57: static double pzero, qzero;
        !            58: static double tpi      = .6366197723675813430755350535e0;
        !            59: static double pio4     = .7853981633974483096156608458e0;
        !            60: static double p1[] = {
        !            61:        0.581199354001606143928050809e21,
        !            62:        -.6672106568924916298020941484e20,
        !            63:        0.2316433580634002297931815435e19,
        !            64:        -.3588817569910106050743641413e17,
        !            65:        0.2908795263834775409737601689e15,
        !            66:        -.1322983480332126453125473247e13,
        !            67:        0.3413234182301700539091292655e10,
        !            68:        -.4695753530642995859767162166e7,
        !            69:        0.2701122710892323414856790990e4,
        !            70: };
        !            71: static double q1[] = {
        !            72:        0.1162398708003212287858529400e22,
        !            73:        0.1185770712190320999837113348e20,
        !            74:        0.6092061398917521746105196863e17,
        !            75:        0.2081661221307607351240184229e15,
        !            76:        0.5243710262167649715406728642e12,
        !            77:        0.1013863514358673989967045588e10,
        !            78:        0.1501793594998585505921097578e7,
        !            79:        0.1606931573481487801970916749e4,
        !            80:        1.0,
        !            81: };
        !            82: static double p2[] = {
        !            83:        -.4435757816794127857114720794e7,
        !            84:        -.9942246505077641195658377899e7,
        !            85:        -.6603373248364939109255245434e7,
        !            86:        -.1523529351181137383255105722e7,
        !            87:        -.1098240554345934672737413139e6,
        !            88:        -.1611616644324610116477412898e4,
        !            89:        0.0,
        !            90: };
        !            91: static double q2[] = {
        !            92:        -.4435757816794127856828016962e7,
        !            93:        -.9934124389934585658967556309e7,
        !            94:        -.6585339479723087072826915069e7,
        !            95:        -.1511809506634160881644546358e7,
        !            96:        -.1072638599110382011903063867e6,
        !            97:        -.1455009440190496182453565068e4,
        !            98:        1.0,
        !            99: };
        !           100: static double p3[] = {
        !           101:        0.3322091340985722351859704442e5,
        !           102:        0.8514516067533570196555001171e5,
        !           103:        0.6617883658127083517939992166e5,
        !           104:        0.1849426287322386679652009819e5,
        !           105:        0.1706375429020768002061283546e4,
        !           106:        0.3526513384663603218592175580e2,
        !           107:        0.0,
        !           108: };
        !           109: static double q3[] = {
        !           110:        0.7087128194102874357377502472e6,
        !           111:        0.1819458042243997298924553839e7,
        !           112:        0.1419460669603720892855755253e7,
        !           113:        0.4002944358226697511708610813e6,
        !           114:        0.3789022974577220264142952256e5,
        !           115:        0.8638367769604990967475517183e3,
        !           116:        1.0,
        !           117: };
        !           118: static double p4[] = {
        !           119:        -.9963753424306922225996744354e23,
        !           120:        0.2655473831434854326894248968e23,
        !           121:        -.1212297555414509577913561535e22,
        !           122:        0.2193107339917797592111427556e20,
        !           123:        -.1965887462722140658820322248e18,
        !           124:        0.9569930239921683481121552788e15,
        !           125:        -.2580681702194450950541426399e13,
        !           126:        0.3639488548124002058278999428e10,
        !           127:        -.2108847540133123652824139923e7,
        !           128:        0.0,
        !           129: };
        !           130: static double q4[] = {
        !           131:        0.5082067366941243245314424152e24,
        !           132:        0.5435310377188854170800653097e22,
        !           133:        0.2954987935897148674290758119e20,
        !           134:        0.1082258259408819552553850180e18,
        !           135:        0.2976632125647276729292742282e15,
        !           136:        0.6465340881265275571961681500e12,
        !           137:        0.1128686837169442121732366891e10,
        !           138:        0.1563282754899580604737366452e7,
        !           139:        0.1612361029677000859332072312e4,
        !           140:        1.0,
        !           141: };
        !           142: 
        !           143: double
        !           144: j1(arg) double arg;{
        !           145:        double xsq, n, d, x;
        !           146:        double sin(), cos(), sqrt();
        !           147:        int i;
        !           148: 
        !           149:        x = arg;
        !           150:        if(x < 0.) x = -x;
        !           151:        if(x > 8.){
        !           152:                asympt(x);
        !           153:                n = x - 3.*pio4;
        !           154:                n = sqrt(tpi/x)*(pzero*cos(n) - qzero*sin(n));
        !           155:                if(arg <0.) n = -n;
        !           156:                return(n);
        !           157:        }
        !           158:        xsq = x*x;
        !           159:        for(n=0,d=0,i=8;i>=0;i--){
        !           160:                n = n*xsq + p1[i];
        !           161:                d = d*xsq + q1[i];
        !           162:        }
        !           163:        return(arg*n/d);
        !           164: }
        !           165: 
        !           166: double
        !           167: y1(arg) double arg;{
        !           168:        double xsq, n, d, x;
        !           169:        double sin(), cos(), sqrt(), log(), j1();
        !           170:        int i;
        !           171: 
        !           172:        x = arg;
        !           173:        if(x <= 0.){
        !           174: #if defined(vax)||defined(tahoe)
        !           175:                extern double infnan();
        !           176:                return(infnan(EDOM));           /* NaN */
        !           177: #else  /* defined(vax)||defined(tahoe) */
        !           178:                return(zero/zero);      /* IEEE machines: invalid operation */
        !           179: #endif /* defined(vax)||defined(tahoe) */
        !           180:        }
        !           181:        if(x > 8.){
        !           182:                asympt(x);
        !           183:                n = x - 3*pio4;
        !           184:                return(sqrt(tpi/x)*(pzero*sin(n) + qzero*cos(n)));
        !           185:        }
        !           186:        xsq = x*x;
        !           187:        for(n=0,d=0,i=9;i>=0;i--){
        !           188:                n = n*xsq + p4[i];
        !           189:                d = d*xsq + q4[i];
        !           190:        }
        !           191:        return(x*n/d + tpi*(j1(x)*log(x)-1./x));
        !           192: }
        !           193: 
        !           194: static
        !           195: asympt(arg) double arg;{
        !           196:        double zsq, n, d;
        !           197:        int i;
        !           198:        zsq = 64./(arg*arg);
        !           199:        for(n=0,d=0,i=6;i>=0;i--){
        !           200:                n = n*zsq + p2[i];
        !           201:                d = d*zsq + q2[i];
        !           202:        }
        !           203:        pzero = n/d;
        !           204:        for(n=0,d=0,i=6;i>=0;i--){
        !           205:                n = n*zsq + p3[i];
        !           206:                d = d*zsq + q3[i];
        !           207:        }
        !           208:        qzero = (8./arg)*(n/d);
        !           209: }

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