The mathematical Gamma function appears in many areas of classical statistics and numerical programming. There are many ways to think about what the Gamma function is. I usually think about Gamma(n) as a function that equals Factorial(n-1) where n can be a real value rather than just an integer.

Here’s part of the Wikipedia entry on the Gamma function:

Implementing Gamma(n) is very difficult. I usually use what is called the Lanczos g=5, n=7 approximation. See the demo code below.

There’s no great moral to this blog post other than that the more techniques you have in your personal toolkit, the more versatile you are when it comes to solving problems.

Mathematics, computer science, and chess all have intrinsic beauty (to me anyway). As I was writing this blog post, the 2021 World Chess Championship was being held in Dubai. The contestants are current world champion Magnus Carlsen (Norway) and the challenger Ian Nepomniachtchi (Russia).

using System;
namespace GammaCSharp
{
  class GammaCSharpProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin Gamma function demo \n");
      Console.WriteLine("For integers, Gamma(n) = (n-1)! \n");
      Console.WriteLine("Example, Gamma(4) = 3! = 3*2*1 = 6 \n");

      int[] nVals = new int[] { 0, 1, 2, 3, 4 };
      for (int i = 0; i < nVals.Length; ++i)
      {
        int n = nVals[i];
        int f = Factorial(n);
        Console.WriteLine("n = " + n.ToString() +
          "     Factorial(n) = " +
          Factorial(n).ToString().PadLeft(5) );
      }

      Console.WriteLine("");

      double[] zVals = new double[] { 0.5, 1.0, 1.5, 2.0, 2.5,
        3.0, 3.5, 4.0, 4.5, 5.0 };
      for (int i = 0; i < zVals.Length; ++i)
      {
        double z = zVals[i];
        double g = Gamma(z);
        Console.WriteLine("z = " + z.ToString("F1") + 
          "    Gamma(z) = " + g.ToString("F4").PadLeft(8));
      }

      Console.WriteLine("\nEnd demo ");
      Console.ReadLine();
    }

    static int Factorial(int n)
    {
      int result = 1;
      for (int i = 1; i <= n; ++i)
        result *= i;
      return result;
    }

    static double Gamma(double z)
    {
      return Math.Exp(LogGamma(z));
    }

    static double LogGamma(double z)
    {
      // Lanczos approximation g=5, n=7
      double[] coef = new double[7] { 1.000000000190015,
         76.18009172947146, -86.50532032941677,
        24.01409824083091, -1.231739572450155,
        0.1208650973866179e-2, -0.5395239384953e-5 };

      double LogSqrtTwoPi = 0.91893853320467274178;

      if (z < 0.5) // Gamma(z) = Pi / (Sin(Pi*z))* Gamma(1-z))
        return Math.Log(Math.PI / Math.Sin(Math.PI * z)) -
          LogGamma(1.0 - z);

      double zz = z - 1.0;
      double b = zz + 5.5; // g + 0.5
      double sum = coef[0];

      for (int i = 1; i < coef.Length; ++i)
        sum += coef[i] / (zz + i);

      return (LogSqrtTwoPi + Math.Log(sum) - b) +
        (Math.Log(b) * (zz + 0.5));
    }
  }
}