Gamma function: implementation?

[John Woods]

In article <430008 at hpqtdla.HP.COM>, bww at hpqtdla.HP.COM (Brian Woodroffe) writes:
> As second question; what is a gamma function anyway?

To answer this first: it is defined as

	GAMMA(z) = integral from 0 to infinity of t^(z-1) * e^(-t) dt

	Real part of z is greater than 0.

It is a generalization of the factorial function (GAMMA(x+1) == x!).  I don't
know what it is good for.

> I am trying to implement a port of a program that calls for a
> "gamma" function.

Find the book "Computer Approximations", by Hart, Cheney, Lawson, Maehly,
Mesztenyi, Rice, Thacher, and Witzgall; ISBN 0-88275-642-7; publisher Robert
E. Krieger Publishing Company (Malabar, Florida) [originally published by
John Wiley & Sons].  "Numerical Recipes in C" is probably also a good book
to check out, but it is at home, so I can't provide a detailed reference.

> for 24-bit significand and a 7 bit exponent; base 2; ie single precision.

For your amusement, I offer the following subroutine, based on the description
found in Hart et al.  I make no promises that it is efficient, especially
accurate in odd cases, or otherwise to the liking of a real, live approximation
weenie.  It works to 6 figure accuracy on the 4 values tested.  If you just
need a quick hack, use it; if you really care about accuracy, get the books
indicated, study the appropriate sections, and code it up yourself (or bash
mine until an approximation weenie would approve).  (In particular, if your
compiler does float arithmetic exclusively with floats, more care might or
might not be needed to handle rounding).

IMPORTANT NOTE:  This routine calculates GAMMA(x).  The UNIX gamma() function
actually calculates ln(GAMMA(x)), which is usually more useful (because
gamma() gets big real fast).  If your application needs the UNIX "gamma()"
function instead of the GAMMA() function, well, check out either reference
(I see a note from Doug Gwyn saying that Numerical Recipes DOES have it).

#include <math.h>	/* for MAXFLOAT */

/*
 * This quick hack calculates gamma for positive reals.  Gamma does have
 * values for negative integers, but checking for integerness is just harder
 * than I am willing to do right now.
 */

float P[] = {	/* Table 5229 */
	 33.0365750763014,
	 16.4294032873095,
	 8.4825775279335,
	 2.1025754321003,
	 .5395052506362,
	 .07829734119881
}, Q[] = {
	 33.0365750692944,
	 2.4620579946557,
	-6.1641661029932,
	 1.0
};

/* Calculate a polynomial given the coefficients from the table. */
/* This does not worry about roundoffs, overflows, or efficiency. */
float Poly(x, coeff, n)
float x;
float coeff[];
int n;
{
	double value;
	int i;
	value = coeff[n];
	for (i = n-1; i >= 0; i--) {
		value = coeff[i] + value * x;
	}
	return value;
}

float fgamma(x) double x;
{
	float mult = 1.0;
	float denom;

	if (x <= 0.0) return MAXFLOAT;

	/* reduce range via recurrence relations *.
	/* want to normalize x to between 2 and 3 */
	if (x < 2.0) {
		/* GAMMA(x) = x GAMMA(x) / x
			    = GAMMA(1+x) / x
		 */
		if (x < 1.0) {
			mult = 1.0 / x;
			x += 1.0;
			/* I am doing this in two steps to make the recurrence
			 * relation more obvious.  A production routine would
			 * do the calculation all at once rather than falling
			 * through
			 */
		}
		mult = mult/x;
		x += 1.0;
	}
	else while (x > 3) {
		mult *= x-1;	/* GAMMA(z+1) = z GAMMA(z) */
		x -= 1.0;
	}
	x -= 2.0;	/* calculate GAMMA(2+x) ~= P(x) / Q(x) */
	denom = Poly(x, P, 5);
	return mult * denom / Poly(x, Q, 3);
}

#ifdef TEST
main() {
	printf("gamma(1/2) is %f\n", fgamma(.5));
	printf("gamma(.2) is %f\n", fgamma(.2));
	printf("gamma(4) is %f\n", fgamma(4.0));
	printf("gamma(5) is %f\n", fgamma(5.0));
}
#endif
-- 
John Woods, Charles River Data Systems, Framingham MA 508-626-1101
...!decvax!frog!john, john at frog.UUCP, ...!mit-eddie!jfw, jfw at eddie.mit.edu