2.11BSD#431 - /usr/src/usr.lib/libom/gamma.c

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

Defined functions

Defined variables

Defined macros