Return to j1.c CVS log

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

Annotation of 43BSDReno/lib/libm/common_source/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.3 (Berkeley) 9/22/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 "mathimpl.h"
        !            52: 
        !            53: #if defined(vax)||defined(tahoe)
        !            54: #include <errno.h>
        !            55: #else  /* defined(vax)||defined(tahoe) */
        !            56: static const double zero = 0.e0;
        !            57: #endif /* defined(vax)||defined(tahoe) */
        !            58: 
        !            59: static double pzero, qzero;
        !            60: 
        !            61: static const double tpi        = .6366197723675813430755350535e0;
        !            62: static const double pio4       = .7853981633974483096156608458e0;
        !            63: static const double p1[] = {
        !            64:        0.581199354001606143928050809e21,
        !            65:        -.6672106568924916298020941484e20,
        !            66:        0.2316433580634002297931815435e19,
        !            67:        -.3588817569910106050743641413e17,
        !            68:        0.2908795263834775409737601689e15,
        !            69:        -.1322983480332126453125473247e13,
        !            70:        0.3413234182301700539091292655e10,
        !            71:        -.4695753530642995859767162166e7,
        !            72:        0.2701122710892323414856790990e4,
        !            73: };
        !            74: static const double q1[] = {
        !            75:        0.1162398708003212287858529400e22,
        !            76:        0.1185770712190320999837113348e20,
        !            77:        0.6092061398917521746105196863e17,
        !            78:        0.2081661221307607351240184229e15,
        !            79:        0.5243710262167649715406728642e12,
        !            80:        0.1013863514358673989967045588e10,
        !            81:        0.1501793594998585505921097578e7,
        !            82:        0.1606931573481487801970916749e4,
        !            83:        1.0,
        !            84: };
        !            85: static const double p2[] = {
        !            86:        -.4435757816794127857114720794e7,
        !            87:        -.9942246505077641195658377899e7,
        !            88:        -.6603373248364939109255245434e7,
        !            89:        -.1523529351181137383255105722e7,
        !            90:        -.1098240554345934672737413139e6,
        !            91:        -.1611616644324610116477412898e4,
        !            92:        0.0,
        !            93: };
        !            94: static const double q2[] = {
        !            95:        -.4435757816794127856828016962e7,
        !            96:        -.9934124389934585658967556309e7,
        !            97:        -.6585339479723087072826915069e7,
        !            98:        -.1511809506634160881644546358e7,
        !            99:        -.1072638599110382011903063867e6,
        !           100:        -.1455009440190496182453565068e4,
        !           101:        1.0,
        !           102: };
        !           103: static const double p3[] = {
        !           104:        0.3322091340985722351859704442e5,
        !           105:        0.8514516067533570196555001171e5,
        !           106:        0.6617883658127083517939992166e5,
        !           107:        0.1849426287322386679652009819e5,
        !           108:        0.1706375429020768002061283546e4,
        !           109:        0.3526513384663603218592175580e2,
        !           110:        0.0,
        !           111: };
        !           112: static const double q3[] = {
        !           113:        0.7087128194102874357377502472e6,
        !           114:        0.1819458042243997298924553839e7,
        !           115:        0.1419460669603720892855755253e7,
        !           116:        0.4002944358226697511708610813e6,
        !           117:        0.3789022974577220264142952256e5,
        !           118:        0.8638367769604990967475517183e3,
        !           119:        1.0,
        !           120: };
        !           121: static const double p4[] = {
        !           122:        -.9963753424306922225996744354e23,
        !           123:        0.2655473831434854326894248968e23,
        !           124:        -.1212297555414509577913561535e22,
        !           125:        0.2193107339917797592111427556e20,
        !           126:        -.1965887462722140658820322248e18,
        !           127:        0.9569930239921683481121552788e15,
        !           128:        -.2580681702194450950541426399e13,
        !           129:        0.3639488548124002058278999428e10,
        !           130:        -.2108847540133123652824139923e7,
        !           131:        0.0,
        !           132: };
        !           133: static const double q4[] = {
        !           134:        0.5082067366941243245314424152e24,
        !           135:        0.5435310377188854170800653097e22,
        !           136:        0.2954987935897148674290758119e20,
        !           137:        0.1082258259408819552553850180e18,
        !           138:        0.2976632125647276729292742282e15,
        !           139:        0.6465340881265275571961681500e12,
        !           140:        0.1128686837169442121732366891e10,
        !           141:        0.1563282754899580604737366452e7,
        !           142:        0.1612361029677000859332072312e4,
        !           143:        1.0,
        !           144: };
        !           145: 
        !           146: static void asympt();
        !           147: 
        !           148: double
        !           149: j1(arg) double arg;{
        !           150:        double xsq, n, d, x;
        !           151:        int i;
        !           152: 
        !           153:        x = arg;
        !           154:        if(x < 0.) x = -x;
        !           155:        if(x > 8.){
        !           156:                asympt(x);
        !           157:                n = x - 3.*pio4;
        !           158:                n = sqrt(tpi/x)*(pzero*cos(n) - qzero*sin(n));
        !           159:                if(arg <0.) n = -n;
        !           160:                return(n);
        !           161:        }
        !           162:        xsq = x*x;
        !           163:        for(n=0,d=0,i=8;i>=0;i--){
        !           164:                n = n*xsq + p1[i];
        !           165:                d = d*xsq + q1[i];
        !           166:        }
        !           167:        return(arg*n/d);
        !           168: }
        !           169: 
        !           170: double
        !           171: y1(arg) double arg;{
        !           172:        double xsq, n, d, x;
        !           173:        int i;
        !           174: 
        !           175:        x = arg;
        !           176:        if(x <= 0.){
        !           177: #if defined(vax)||defined(tahoe)
        !           178:                return(infnan(EDOM));           /* NaN */
        !           179: #else  /* defined(vax)||defined(tahoe) */
        !           180:                return(zero/zero);      /* IEEE machines: invalid operation */
        !           181: #endif /* defined(vax)||defined(tahoe) */
        !           182:        }
        !           183:        if(x > 8.){
        !           184:                asympt(x);
        !           185:                n = x - 3*pio4;
        !           186:                return(sqrt(tpi/x)*(pzero*sin(n) + qzero*cos(n)));
        !           187:        }
        !           188:        xsq = x*x;
        !           189:        for(n=0,d=0,i=9;i>=0;i--){
        !           190:                n = n*xsq + p4[i];
        !           191:                d = d*xsq + q4[i];
        !           192:        }
        !           193:        return(x*n/d + tpi*(j1(x)*log(x)-1./x));
        !           194: }
        !           195: 
        !           196: static void
        !           197: asympt(arg) double arg;{
        !           198:        double zsq, n, d;
        !           199:        int i;
        !           200:        zsq = 64./(arg*arg);
        !           201:        for(n=0,d=0,i=6;i>=0;i--){
        !           202:                n = n*zsq + p2[i];
        !           203:                d = d*zsq + q2[i];
        !           204:        }
        !           205:        pzero = n/d;
        !           206:        for(n=0,d=0,i=6;i>=0;i--){
        !           207:                n = n*zsq + p3[i];
        !           208:                d = d*zsq + q3[i];
        !           209:        }
        !           210:        qzero = (8./arg)*(n/d);
        !           211: }

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