Return to besj1.c CVS log

Up to [Research Unix] / researchv10dc / libc / math

Annotation of researchv10dc/libc/math/besj1.c, revision 1.1

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

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