Return to j0.c CVS log

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

Annotation of 43BSDReno/lib/libm/common_source/j0.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[] = "@(#)j0.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 zero
        !            15: 
        !            16:        j0(x) returns the value of J0(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 J0 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:        #5849 (19.22D)
        !            31:        #6549 (19.25D)
        !            32:        #6949 (19.41D)
        !            33: 
        !            34:        y0(x) returns the value of Y0(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, j0.
        !            41: 
        !            42:        The values of Y0 have not been checked
        !            43:        to more than ten places.
        !            44: 
        !            45:        Coefficients are from Hart & Cheney.
        !            46:        #6245 (18.78D)
        !            47:        #6549 (19.25D)
        !            48:        #6949 (19.41D)
        !            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.4933787251794133561816813446e21,
        !            65:        -.1179157629107610536038440800e21,
        !            66:        0.6382059341072356562289432465e19,
        !            67:        -.1367620353088171386865416609e18,
        !            68:        0.1434354939140344111664316553e16,
        !            69:        -.8085222034853793871199468171e13,
        !            70:        0.2507158285536881945555156435e11,
        !            71:        -.4050412371833132706360663322e8,
        !            72:        0.2685786856980014981415848441e5,
        !            73: };
        !            74: static const double q1[] = {
        !            75:        0.4933787251794133562113278438e21,
        !            76:        0.5428918384092285160200195092e19,
        !            77:        0.3024635616709462698627330784e17,
        !            78:        0.1127756739679798507056031594e15,
        !            79:        0.3123043114941213172572469442e12,
        !            80:        0.6699987672982239671814028660e9,
        !            81:        0.1114636098462985378182402543e7,
        !            82:        0.1363063652328970604442810507e4,
        !            83:        1.0
        !            84: };
        !            85: static const double p2[] = {
        !            86:        0.5393485083869438325262122897e7,
        !            87:        0.1233238476817638145232406055e8,
        !            88:        0.8413041456550439208464315611e7,
        !            89:        0.2016135283049983642487182349e7,
        !            90:        0.1539826532623911470917825993e6,
        !            91:        0.2485271928957404011288128951e4,
        !            92:        0.0,
        !            93: };
        !            94: static const double q2[] = {
        !            95:        0.5393485083869438325560444960e7,
        !            96:        0.1233831022786324960844856182e8,
        !            97:        0.8426449050629797331554404810e7,
        !            98:        0.2025066801570134013891035236e7,
        !            99:        0.1560017276940030940592769933e6,
        !           100:        0.2615700736920839685159081813e4,
        !           101:        1.0,
        !           102: };
        !           103: static const double p3[] = {
        !           104:        -.3984617357595222463506790588e4,
        !           105:        -.1038141698748464093880530341e5,
        !           106:        -.8239066313485606568803548860e4,
        !           107:        -.2365956170779108192723612816e4,
        !           108:        -.2262630641933704113967255053e3,
        !           109:        -.4887199395841261531199129300e1,
        !           110:        0.0,
        !           111: };
        !           112: static const double q3[] = {
        !           113:        0.2550155108860942382983170882e6,
        !           114:        0.6667454239319826986004038103e6,
        !           115:        0.5332913634216897168722255057e6,
        !           116:        0.1560213206679291652539287109e6,
        !           117:        0.1570489191515395519392882766e5,
        !           118:        0.4087714673983499223402830260e3,
        !           119:        1.0,
        !           120: };
        !           121: static const double p4[] = {
        !           122:        -.2750286678629109583701933175e20,
        !           123:        0.6587473275719554925999402049e20,
        !           124:        -.5247065581112764941297350814e19,
        !           125:        0.1375624316399344078571335453e18,
        !           126:        -.1648605817185729473122082537e16,
        !           127:        0.1025520859686394284509167421e14,
        !           128:        -.3436371222979040378171030138e11,
        !           129:        0.5915213465686889654273830069e8,
        !           130:        -.4137035497933148554125235152e5,
        !           131: };
        !           132: static const double q4[] = {
        !           133:        0.3726458838986165881989980e21,
        !           134:        0.4192417043410839973904769661e19,
        !           135:        0.2392883043499781857439356652e17,
        !           136:        0.9162038034075185262489147968e14,
        !           137:        0.2613065755041081249568482092e12,
        !           138:        0.5795122640700729537480087915e9,
        !           139:        0.1001702641288906265666651753e7,
        !           140:        0.1282452772478993804176329391e4,
        !           141:        1.0,
        !           142: };
        !           143: 
        !           144: static void asympt();
        !           145: 
        !           146: double
        !           147: j0(arg) double arg;{
        !           148:        double argsq, n, d;
        !           149:        int i;
        !           150: 
        !           151:        if(arg < 0.) arg = -arg;
        !           152:        if(arg > 8.){
        !           153:                asympt(arg);
        !           154:                n = arg - pio4;
        !           155:                return(sqrt(tpi/arg)*(pzero*cos(n) - qzero*sin(n)));
        !           156:        }
        !           157:        argsq = arg*arg;
        !           158:        for(n=0,d=0,i=8;i>=0;i--){
        !           159:                n = n*argsq + p1[i];
        !           160:                d = d*argsq + q1[i];
        !           161:        }
        !           162:        return(n/d);
        !           163: }
        !           164: 
        !           165: double
        !           166: y0(arg) double arg;{
        !           167:        double argsq, n, d;
        !           168:        int i;
        !           169: 
        !           170:        if(arg <= 0.){
        !           171: #if defined(vax)||defined(tahoe)
        !           172:                return(infnan(EDOM));           /* NaN */
        !           173: #else  /* defined(vax)||defined(tahoe) */
        !           174:                return(zero/zero);      /* IEEE machines: invalid operation */
        !           175: #endif /* defined(vax)||defined(tahoe) */
        !           176:        }
        !           177:        if(arg > 8.){
        !           178:                asympt(arg);
        !           179:                n = arg - pio4;
        !           180:                return(sqrt(tpi/arg)*(pzero*sin(n) + qzero*cos(n)));
        !           181:        }
        !           182:        argsq = arg*arg;
        !           183:        for(n=0,d=0,i=8;i>=0;i--){
        !           184:                n = n*argsq + p4[i];
        !           185:                d = d*argsq + q4[i];
        !           186:        }
        !           187:        return(n/d + tpi*j0(arg)*log(arg));
        !           188: }
        !           189: 
        !           190: static void
        !           191: asympt(arg) double arg;{
        !           192:        double zsq, n, d;
        !           193:        int i;
        !           194:        zsq = 64./(arg*arg);
        !           195:        for(n=0,d=0,i=6;i>=0;i--){
        !           196:                n = n*zsq + p2[i];
        !           197:                d = d*zsq + q2[i];
        !           198:        }
        !           199:        pzero = n/d;
        !           200:        for(n=0,d=0,i=6;i>=0;i--){
        !           201:                n = n*zsq + p3[i];
        !           202:                d = d*zsq + q3[i];
        !           203:        }
        !           204:        qzero = (8./arg)*(n/d);
        !           205: }

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