Return to lgamma.c CVS log

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

Annotation of 43BSDTahoe/usr.lib/libm/lgamma.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[] = "@(#)lgamma.c   5.2 (Berkeley) 4/29/88";
        !             9: #endif /* not lint */
        !            10: 
        !            11: /*
        !            12:        C program for floating point log Gamma function
        !            13: 
        !            14:        lgamma(x) computes the log of the absolute
        !            15:        value of the Gamma function.
        !            16:        The sign of the Gamma function is returned in the
        !            17:        external quantity signgam.
        !            18: 
        !            19:        The coefficients for expansion around zero
        !            20:        are #5243 from Hart & Cheney; for expansion
        !            21:        around infinity they are #5404.
        !            22: 
        !            23:        Calls log, floor and sin.
        !            24: */
        !            25: 
        !            26: #include <math.h>
        !            27: #if defined(vax)||defined(tahoe)
        !            28: #include <errno.h>
        !            29: #endif /* defined(vax)||defined(tahoe) */
        !            30: int    signgam = 0;
        !            31: static double goobie   = 0.9189385332046727417803297;  /* log(2*pi)/2 */
        !            32: static double pi       = 3.1415926535897932384626434;
        !            33: 
        !            34: #define M 6
        !            35: #define N 8
        !            36: static double p1[] = {
        !            37:        0.83333333333333101837e-1,
        !            38:        -.277777777735865004e-2,
        !            39:        0.793650576493454e-3,
        !            40:        -.5951896861197e-3,
        !            41:        0.83645878922e-3,
        !            42:        -.1633436431e-2,
        !            43: };
        !            44: static double p2[] = {
        !            45:        -.42353689509744089647e5,
        !            46:        -.20886861789269887364e5,
        !            47:        -.87627102978521489560e4,
        !            48:        -.20085274013072791214e4,
        !            49:        -.43933044406002567613e3,
        !            50:        -.50108693752970953015e2,
        !            51:        -.67449507245925289918e1,
        !            52:        0.0,
        !            53: };
        !            54: static double q2[] = {
        !            55:        -.42353689509744090010e5,
        !            56:        -.29803853309256649932e4,
        !            57:        0.99403074150827709015e4,
        !            58:        -.15286072737795220248e4,
        !            59:        -.49902852662143904834e3,
        !            60:        0.18949823415702801641e3,
        !            61:        -.23081551524580124562e2,
        !            62:        0.10000000000000000000e1,
        !            63: };
        !            64: 
        !            65: double
        !            66: lgamma(arg)
        !            67: double arg;
        !            68: {
        !            69:        double log(), pos(), neg(), asym();
        !            70: 
        !            71:        signgam = 1.;
        !            72:        if(arg <= 0.) return(neg(arg));
        !            73:        if(arg > 8.) return(asym(arg));
        !            74:        return(log(pos(arg)));
        !            75: }
        !            76: 
        !            77: static double
        !            78: asym(arg)
        !            79: double arg;
        !            80: {
        !            81:        double log();
        !            82:        double n, argsq;
        !            83:        int i;
        !            84: 
        !            85:        argsq = 1./(arg*arg);
        !            86:        for(n=0,i=M-1; i>=0; i--){
        !            87:                n = n*argsq + p1[i];
        !            88:        }
        !            89:        return((arg-.5)*log(arg) - arg + goobie + n/arg);
        !            90: }
        !            91: 
        !            92: static double
        !            93: neg(arg)
        !            94: double arg;
        !            95: {
        !            96:        double t;
        !            97:        double log(), sin(), floor(), pos();
        !            98: 
        !            99:        arg = -arg;
        !           100:      /*
        !           101:       * to see if arg were a true integer, the old code used the
        !           102:       * mathematically correct observation:
        !           103:       * sin(n*pi) = 0 <=> n is an integer.
        !           104:       * but in finite precision arithmetic, sin(n*PI) will NEVER
        !           105:       * be zero simply because n*PI is a rational number.  hence
        !           106:       *        it failed to work with our newer, more accurate sin()
        !           107:       * which uses true pi to do the argument reduction...
        !           108:       *        temp = sin(pi*arg);
        !           109:       */
        !           110:        t = floor(arg);
        !           111:        if (arg - t  > 0.5e0)
        !           112:            t += 1.e0;                          /* t := integer nearest arg */
        !           113: #if defined(vax)||defined(tahoe)
        !           114:        if (arg == t) {
        !           115:            extern double infnan();
        !           116:            return(infnan(ERANGE));             /* +INF */
        !           117:        }
        !           118: #endif /* defined(vax)||defined(tahoe) */
        !           119:        signgam = (int) (t - 2*floor(t/2));     /* signgam =  1 if t was odd, */
        !           120:                                                /*            0 if t was even */
        !           121:        signgam = signgam - 1 + signgam;        /* signgam =  1 if t was odd, */
        !           122:                                                /*           -1 if t was even */
        !           123:        t = arg - t;                            /*  -0.5 <= t <= 0.5 */
        !           124:        if (t < 0.e0) {
        !           125:            t = -t;
        !           126:            signgam = -signgam;
        !           127:        }
        !           128:        return(-log(arg*pos(arg)*sin(pi*t)/pi));
        !           129: }
        !           130: 
        !           131: static double
        !           132: pos(arg)
        !           133: double arg;
        !           134: {
        !           135:        double n, d, s;
        !           136:        register i;
        !           137: 
        !           138:        if(arg < 2.) return(pos(arg+1.)/arg);
        !           139:        if(arg > 3.) return((arg-1.)*pos(arg-1.));
        !           140: 
        !           141:        s = arg - 2.;
        !           142:        for(n=0,d=0,i=N-1; i>=0; i--){
        !           143:                n = n*s + p2[i];
        !           144:                d = d*s + q2[i];
        !           145:        }
        !           146:        return(n/d);
        !           147: }

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