4.3BSD - /usr/src/usr.lib/libm/lgamma.c

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

Defined functions

Defined variables

Defined macros