Return to j1.c CVS log

Up to [CSRG BSD Unix] / 43BSDTahoe / old / libom

Annotation of 43BSDTahoe/old/libom/j1.c, revision 1.1

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

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