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