Return to j1.c CVS log

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

Annotation of 43BSD/usr.lib/libm/j1.c, revision 1.1

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

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