Return to gamma.c CVS log

Up to [CSRG BSD Unix] / 43BSD / contrib / apl / src

Annotation of 43BSD/contrib/apl/src/gamma.c, revision 1.1.1.1

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

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