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

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.