bww@hpqtdla.HP.COM (Brian Woodroffe) (09/15/89)
I am trying to implement a port of a program that calls for a "gamma" function. Unfortunately the compiler / math's library I am using does not have such a function. Can any one post, or mail me a copy of a PD - C program to implement a gamma function? As second question; what is a gamma function anyway? (I understand that a factorial is the degenerate integer case.) And what is a gamma function used for and where? Thank you, Brian Woodroffe. PS: I am interrested in an implementation that is accurate enough for 24-bit significand and a 7 bit exponent; base 2; ie single precision. +--------------------------------+-----------------------------------+ | Brian Woodroffe | | | Hewlett Packard Ltd | ARPA: bww@hpsqf | | Queensferry Telecomms Division | bww%hpsqf@hplabs.hp.com | | South Queensferry | UUCP: ..!hplabs!hpsqf!bww | | West Lothian | JANET: bww%hpqtdla@hpl.hp.co.uk | | Scotland EH30 9TG. | PHONE: +44-31-331-7234 | +--------------------------------+-----------------------------------+
gwyn@smoke.BRL.MIL (Doug Gwyn) (09/17/89)
In article <430008@hpqtdla.HP.COM> bww@hpqtdla.HP.COM (Brian Woodroffe) writes: >I am trying to implement a port of a program that calls for a >"gamma" function. Unfortunately the compiler / math's library >I am using does not have such a function. Check for a "log of gamma" function. Because gamma grows at a rapid rate, that's generally a more useful function (usually called lgamma()). >As second question; what is a gamma function anyway? (I understand >that a factorial is the degenerate integer case.) And what is a >gamma function used for and where? Heavens, they have lots of applications. Look in a few textbooks. For an implementation, check out Chapter 6 of Press, Flannery, Teukolsky, & Vettering's "Numerical Recipes in C".
mccaugh@s.cs.uiuc.edu (09/21/89)
Brian:
I am posting this in case mailing fails (also in case anyone wants to correct):
#include <math.h>
/* Gamma - xx = argument, */
/* ier = return-code (0: no error; 1: negative; 2: overflow */
double Gamma (xx,ier)
double xx;
int *ier;
{
double x,y, Gx,Gy, err;
#define MAX_FLOAT 1.0E75
if(xx > 57.0) { *ier = 2; /* overflow */
return(MAX_FLOAT);
}
/* else xx <= 57.0 */
x = xx;
err = 1.0E-6; /* margin of error */
*ier = 0;
Gx = 1.0;
if(x > 2.0) { do { x = x - 1.0;
Gx = Gx*x;
}
while(x > 2.0);
/* x <= 2.0 */
y = x - 1.0;
Gy = 1.0 - 0.57710166*(y-1) + 0.98585399*(y-1)^2 -
0.87642182*(y-1)^3 + 0.83282120*(y-1)^4 -
0.56847290*(y-1)^5 + 0.25482049*(y-1)^6 -
0.05149930*(y-1)^7;
Gx = Gx*Gy;
return(Gx);
}
/* else x <= 2.0 */
if(x == 1.0) return(Gx);
if(x > 1.0) { y = x - 1.0;
Gy = 1.0 - 0.57710166*(y-1) + 0.98585399*(y-1)^2 -
0.87642182*(y-1)^3 + 0.83282120*(y-1)^4 -
0.56847290*(y-1)^5 + 0.25482049*(y-1)^6 -
0.05149930*(y-1)^7;
Gx = Gx*Gy;
return(Gx);
}
/* else x < 1.0 */
if(x > err) { do { Gx = Gx/x;
x = x + 1.0;
}
while(x <= 1.0);
/* x > 1.0 */
y = x - 1.0;
Gy = 1.0 - 0.57710166*(y-1) + 0.98585399*(y-1)^2 -
0.87642182*(y-1)^3 + 0.83282120*(y-1)^4 -
0.56847290*(y-1)^5 + 0.25482049*(y-1)^6 -
0.05149930*(y-1)^7;
Gx = Gx*Gy;
return(Gx);
}
/* else x <= err */
y = floor(x) - x;
if(fabs(y) <= err) { *ier = 1;
return(Gx);
}
/* else fabs(y) > err */
if(1-y <= err) { *ier = 1;
return(Gx);
}
/* else 1-y > err */
while(x <= 1.0) { Gx = Gx/x;
x = x + 1.0;
}
/* x > 1.0 */
y = x - 1.0;
Gy = 1.0 - 0.57710166*(y-1) + 0.98585399*(y-1)^2 -
0.87642182*(y-1)^3 + 0.83282120*(y-1)^4 -
0.56847290*(y-1)^5 + 0.25482049*(y-1)^6 -
0.05149930*(y-1)^7;
Gx = Gx*Gy;
return(Gx);
}
My home-brewed C takes exponentiation ('^') while other C's may not; other-
wise, this should work. Hope this helps, and good luck with your porting!
Scott McCaughrin
University of Illinois
Urbana, Illinois.john@frog.UUCP (John Woods) (09/22/89)
In article <430008@hpqtdla.HP.COM>, bww@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@frog.UUCP, ...!mit-eddie!jfw, jfw@eddie.mit.edu