Return to gamma.c CVS log

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

Annotation of researchv9/libc/math/gamma.c, revision 1.1

1.1     ! root        1: /*
        !             2:        C program for floating point log gamma function
        !             3: 
        !             4:        gamma(x) computes the log of the absolute
        !             5:        value of the gamma function.
        !             6:        The sign of the gamma function is returned in the
        !             7:        external quantity signgam.
        !             8: 
        !             9:        The coefficients for expansion around zero
        !            10:        are #5243 from Hart & Cheney; for expansion
        !            11:        around infinity they are #5404.
        !            12: 
        !            13:        Calls log and sin.
        !            14: */
        !            15: 
        !            16: #include <errno.h>
        !            17: #include <math.h>
        !            18: 
        !            19: int    errno;
        !            20: int    signgam = 0;
        !            21: static double goobie   = 0.9189385332046727417803297;
        !            22: static double pi       = 3.1415926535897932384626434;
        !            23: 
        !            24: #define M 6
        !            25: #define N 8
        !            26: static double p1[] = {
        !            27:        0.83333333333333101837e-1,
        !            28:        -.277777777735865004e-2,
        !            29:        0.793650576493454e-3,
        !            30:        -.5951896861197e-3,
        !            31:        0.83645878922e-3,
        !            32:        -.1633436431e-2,
        !            33: };
        !            34: static double p2[] = {
        !            35:        -.42353689509744089647e5,
        !            36:        -.20886861789269887364e5,
        !            37:        -.87627102978521489560e4,
        !            38:        -.20085274013072791214e4,
        !            39:        -.43933044406002567613e3,
        !            40:        -.50108693752970953015e2,
        !            41:        -.67449507245925289918e1,
        !            42:        0.0,
        !            43: };
        !            44: static double q2[] = {
        !            45:        -.42353689509744090010e5,
        !            46:        -.29803853309256649932e4,
        !            47:        0.99403074150827709015e4,
        !            48:        -.15286072737795220248e4,
        !            49:        -.49902852662143904834e3,
        !            50:        0.18949823415702801641e3,
        !            51:        -.23081551524580124562e2,
        !            52:        0.10000000000000000000e1,
        !            53: };
        !            54: 
        !            55: double
        !            56: gamma(arg)
        !            57: double arg;
        !            58: {
        !            59:        double log(), pos(), neg(), asym();
        !            60: 
        !            61:        signgam = 1.;
        !            62:        if(arg <= 0.) return(neg(arg));
        !            63:        if(arg > 8.) return(asym(arg));
        !            64:        return(log(pos(arg)));
        !            65: }
        !            66: 
        !            67: static double
        !            68: asym(arg)
        !            69: double arg;
        !            70: {
        !            71:        double log();
        !            72:        double n, argsq;
        !            73:        int i;
        !            74: 
        !            75:        argsq = 1./(arg*arg);
        !            76:        for(n=0,i=M-1; i>=0; i--){
        !            77:                n = n*argsq + p1[i];
        !            78:        }
        !            79:        return((arg-.5)*log(arg) - arg + goobie + n/arg);
        !            80: }
        !            81: 
        !            82: static double
        !            83: neg(arg)
        !            84: double arg;
        !            85: {
        !            86:        double temp;
        !            87:        double log(), sin(), pos();
        !            88: 
        !            89:        arg = -arg;
        !            90:        temp = sin(pi*arg);
        !            91:        if(temp == 0.) {
        !            92:                errno = EDOM;
        !            93:                return(HUGE);
        !            94:        }
        !            95:        if(temp < 0.) temp = -temp;
        !            96:        else signgam = -1;
        !            97:        return(-log(arg*pos(arg)*temp/pi));
        !            98: }
        !            99: 
        !           100: static double
        !           101: pos(arg)
        !           102: double arg;
        !           103: {
        !           104:        double n, d, s;
        !           105:        register i;
        !           106: 
        !           107:        if(arg < 2.) return(pos(arg+1.)/arg);
        !           108:        if(arg > 3.) return((arg-1.)*pos(arg-1.));
        !           109: 
        !           110:        s = arg - 2.;
        !           111:        for(n=0,d=0,i=N-1; i>=0; i--){
        !           112:                n = n*s + p2[i];
        !           113:                d = d*s + q2[i];
        !           114:        }
        !           115:        return(n/d);
        !           116: }

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