## The Log Gamma Function with C#

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;
}
```