I was doing some numerical programming recently and needed a C# implementation of the incomplete gamma function. This in turn required a C# implementation of the log-gamma function. (Note: the gamma function is a cousin to the factorial function. The incomplete gamma function is a cousin to the gamma function. The incomplete gamma function can be used to compute probabilities related to the chi-square statistic, which is used when comparing a set of observed values with a set of expected values.)

So, I pulled up the online version of the famous “Handbook of Mathematical Functions” by Abramowitz and Stegun (often just called A&S) and starting coding. I wrote three versions of a LogGamma function. The first, and most accurate, uses an amazing technique called the Lanczos algorithm. There are several variations of the Lanczos algorithm and I used the one called “g=5, n=6/7”. The second implementation of LogGamma used a continuing fraction equation (A&S 6.1.48). The third implementation used a series expansion equation (A&S 6.1.41).

The accuracy of the methods depends to some extent on the value of the input, x. For small values of x (less than 2), the Lanczos algorithm implementation was accurate to 10 decimals. The continuing fraction version was accurate to about 4 decimals. The series expansion version was pretty bad and only accurate to about 3 decimals (but was much better for larger values of x). In the screenshot below, I compare computed values against valuews from a table in A&S.

There’s no moral to this blog post. But if you ever need a C# implementation of the log-gamma function, you can use the code below.

static double LogGammaLanczos(double x) { // Log of Gamma from Lanczos with g=5, n=6/7 // not in A & S double[] coef = new double[6] { 76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 0.1208650973866179E-2, -0.5395239384953E-5 }; double LogSqrtTwoPi = 0.91893853320467274178; double denom = x + 1; double y = x + 5.5; double series = 1.000000000190015; for (int i = 0; i < 6; ++i) { series += coef[i] / denom; denom += 1.0; } return (LogSqrtTwoPi + (x + 0.5) * Math.Log(y) - y + Math.Log(series / x)); }

static double LogGammaContinued(double x) { // A & S eq. 6.1.48 (continuing fraction) double a0 = 1.0 / 12; double a1 = 1.0 / 30; double a2 = 53.0 / 210; double a3 = 195.0 / 371; double a4 = 22999.0 / 22737; double a5 = 29944523.0 / 19733142; double a6 = 109535241009.0 / 48264275462; double t6 = a6 / x; double t5 = a5 / (x + t6); double t4 = a4 / (x + t5); double t3 = a3 / (x + t4); double t2 = a2 / (x + t3); double t1 = a1 / (x + t2); double t0 = a0 / (x + t1); double result = t0 - x + ((x - 0.5) * Math.Log(x)) + (0.5 * Math.Log(2 * Math.PI)); return result; }

static double LogGammaSeries(double z) { // A & S 6.1.41 (Stirling's approximation) double x1 = (z - 0.5) * Math.Log(z); double x3 = 0.5 * Math.Log(2 * Math.PI); double x4 = 1 / (12 * z); double x5 = 1 / (360 * z * z * z); double x6 = 1 / (1260 * z * z * z * z * z); double x7 = 1 / (1680 * z * z * z * z * z * z * z); // more terms possible return x1 - z + x3 + x4 - x5 + x6 - x7; }