2.11BSD#431 - /usr/src/new/PORT/apl/src/gamma.c

   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: }

Defined functions

Defined variables

Defined macros