Let’s flex our mathematical legs a bit by implementing an approximation to the gamma function Γ(z). In this article, I’ll introduce what the gamma function is, what it’s for, and how to actually calculate it, using a well known approximation method.
This topic leans more heavily into mathematics than usual for this blog. Some of the concepts may be new to you, and I regret I cannot give them a proper treatment in a single article. I’ll do my best to give you the flavor of how everything works, and I have provided some resources at the end should you wish to learn more.
What’s the gamma function? What’s it for?
You’re probably familiar with integer factorials already. Here’s a simple definition:
\(n! = 1 \times 2 \times 3 \times \cdots \times (n-2) \times (n-1) \times n\)
For example, \(5! = 1 \times 2 \times 3 \times 4 \times 5 = 120 \). However, that only works for positive integers. The gamma function lets us work with most positive and negative real numbers, and even complex numbers. For positive integers, the gamma function is very closely related to ordinary factorials:
\(\Gamma(z) = (z – 1)!\)
Why is it (z – 1)! instead of just z! ? As H. M. Edwards puts it—in a footnote, no less:
Unfortunately, Legendre subsequently introduced the notation of \(\Gamma(s) \) for \(\Pi(s-1). \) Legendre’s reasons for considering (n – 1)! instead of n! are obscure (perhaps he felt it was more natural to have the first pole occur at s = 0 rather than at s = -1), but, whatever the reason, this notation prevailed in France and, by the end of the nineteenth century, in the rest of the world as well. Gauss’s original notation appears to me to be much more natural … .
H. M. Edwards in Riemann’s Zeta Function
A couple of centuries later, it seems we’re stuck with \(\Gamma(z) = (z – 1)! \). The gamma function is defined by the following integral:
\({\displaystyle \Gamma(z)=\int _{0}^{\infty }t^{z-1}e^{-t}\,dt}\)
This is one of two integrals to bear the name “Euler’s integral.” Hopefully that gives you some additional sense of just how important it is.
Complex numbers, far too quickly:
Recall that complex numbers can be represented as follows:
\(z = a + bi\)
Where a and b are the real and imaginary parts to the number. i is the imaginary unit, satisfying i2 = −1, though for our purposes today we can think of it as a suffix that means “this is the imaginary part”. For example, the number 3 – 5i is a complex number with a real part of 3 and an imaginary part of -5.
A real number is a complex number with an imaginary part of zero. For example, 3.2 could be written 3.2 + 0i.
The Lanczos approximation
The above integral is beautiful, but is not practical if we want to compute a decimal approximation of a gamma value. To that end, we’ll look at the Lanczos method, developed by the Hungarian-Irish mathematician and physicist Cornelius Lanczos after World War II.
Lanczos served as an assistant to Einstein between 1928 and 1929, having earned some praise for his Ph.D. thesis before that.
Below is the formula most typically given for the Lanczos approximation, which is valid for \(\Re(z) > \frac{1}{2} \), that is, numbers with a real part greater than 0.5:
\({\displaystyle \Gamma (z+1)=
\frac{{\sqrt {2\pi }}{\left(z+g+{\frac {1}{2}}\right)}^{z+\frac{1}{2}}}
{e^{z+g+\frac{1}{2}}}
A_{g}(z)}
\)
With \(A_g(z) \) defined as:
\({\displaystyle A_{g}(z)={\frac {1}{2}}p_{0}(g)+p_{1}(g){\frac {z}{z+1}}+p_{2}(g){\frac {z(z-1)}{(z+1)(z+2)}}+\cdots}\)
The constant g is an arbitrary value that one can tweak to get a faster or slower convergence. I’ve gone with a value of 7, which is common. The \(p_n(g) \) values are coefficients generated through a process I will not cover here. There are other published works that provide these coefficients for us, so I’ve opted to use those. Particular thanks to Robert Munafo’s page listing several coefficient tables:
# Credit: https://mrob.com/pub/ries/lanczos-gamma.html
my @p = (
0.99999999999980993227684700473478,
676.520368121885098567009190444019,
-1259.13921672240287047156078755283,
771.3234287776530788486528258894,
-176.61502916214059906584551354,
12.507343278686904814458936853,
-0.13857109526572011689554707,
9.984369578019570859563e-6,
1.50563273514931155834e-7,
);
Optimizing for accuracy
Implementing the above formulas directly can lead to roundoff errors with floating point arithmetic, so the Lanczos approximation is very often written in a partial fraction form that also happens to be a little easier to compute:
\({\displaystyle \Gamma(z) = \sqrt{2 \pi} \frac{(z + g – \frac{1}{2})^{z-\frac{1}{2}}}{e^{z+g-\frac{1}{2}}} A_g(z)}, \text{or} \\
\Gamma(z) = \sqrt{2 \pi} (z + g – \frac{1}{2})^{z-\frac{1}{2}} e^{-(z+g-\frac{1}{2})} A_g(z)\)
With \(A_g(z) \) defined as the much simpler sum:
\({\displaystyle A_{g}(z)= p_{0}(g) + \sum_{i=1}^{n} { \frac{p_{i}(g)}{z + l} } }\)
Other Approximations
Other approximations to the gamma function exist, which will appear similar to Lanczos’s method. Some authors even refer to them incorrectly as a Lanczos approximation. Of these, the Spouge approximation and Stirling approximation are probably the two most common.
Reflection Formula
Remember, the above formula is only valid for \(\Re(z) > \frac{1}{2} \). That is, if the real part of z is less than or equal to one half, it won’t work. Fortunately, for \(\Re(z) \leq \frac{1}{2} \), we can use the following reflection formula:
\({\displaystyle \Gamma(z) \cdot \Gamma(1-z) = \frac{\pi}{\sin \pi z} } \quad \text{for } \Re(z) \leq \frac{1}{2}\)
This can be rearranged to:
\({\displaystyle \Gamma(z) = \frac{\pi}{\sin \pi z \times \Gamma(1-z) } } \quad \text{for } \Re(z) \leq \frac{1}{2}\)
We’re finally there. We finally have a form that we can use to get a discrete answer, and after all that math, here is the result, which ends up only needing a few lines of code:
# Lanczos approximation of the gamma function.
# $z may be real or complex (Math::Complex)
sub gamma($z) {
return pi / (sin($z * pi) * gamma(1 - $z)) if Re($z) < 0.5;
my $A = $p[0] + sum map { $p[$_] / ($z + $_ - 1) } 1..$#p;
my $t = $z + g - 0.5;
sqrt(2*pi) * $A * $t**($z-0.5) * exp(-$t);
}
In the code above, we start with the reflection formula (which may recurse). $A is the \(A_g(z) \) sum, $t is convenience for the reused $z + g - 0.5 part, and the last line ties the rest of the main formula together.
That’s it! Thanks to the overloaded operators in Math::Complex, our gamma($z) will work with real or complex numbers without doing anything special.
ch-2.pl full source.
Unit Testing
I only include unit tests in these blog posts when I feel they offer some value to the reader. In this case, I can’t just use is(), because I’m comparing floating point values which may vary slightly. As such, I need to ensure the gamma values match expected values within a small tolerance, which I’ve set as a constant (tol => 1e-7).
I’m writing my own tolerance check to avoid a non-core dependency. Normally I would use Test::Number::Delta within => 1e-7.
Then I could write delta_ok(Re(gamma($z)), Re($exp)).
I needed to test both real and complex inputs/outputs, so I defined a sub that will handle both:
sub gamma_ok($z,$exp) {
my $gamma = gamma(cplx($z));
ok(( abs(Re($gamma) - Re($exp)) <= tol
and abs(Im($gamma) - Im($exp)) <= tol ),
"gamma($z) = $exp within " .tol)
or diag " got $gamma";
}
Here, $z is the input (which can be a real or complex number, or complex string accepted by the Math::Complex::cplx() function, such as '5-3i'. $exp is the expected value, which can also be real or complex. I will now define my desired tests in a hash and loop over them:
use Test::More;
my %tests = (
# Real tests
1 => 1,
2 => 1,
3 => 2,
4 => 6,
5 => 24,
0.5 => sqrt(pi),
1.5 => 0.5 * sqrt(pi),
-0.5 => -2 * sqrt(pi),
# Complex tests
'1-1i' => 0.4980156681 + 0.1549498283*i,
'0.5+0.5i' => 0.8181639995 - 0.7633138287*i,
'5+3i' => 0.0160418827 - 9.4332932898*i,
'5-3i' => 0.0160418827 + 9.4332932898*i,
);
gamma_ok($_, $tests{$_}) for sort keys %tests;
done_testing;
A note on negative z values
Remember, the Lanczos function is a representation of the gamma function for z > 0. However, the function exists for all z except negative integers and zero, where it is undefined. Here is a graph I very roughly plotted with the SVG code from a few weeks ago.
Near those undefined (asymptotic) points, my gamma($z) function will return very large or small values depending on whether the limit approaches positive or negative infinity from that direction.
You might think I could just return Inf() if $z is a negative integer or zero, but there is no good way to tell if someone really wants to know the value of gamma(-2), or whether they meant gamma(-1.999999999999999999999) and Perl rounded up (or vice versa), at least not without making a lot of potentially surprising assumptions about the input. Therefore, I believe it makes the most sense to trust the caller to know what they want, and handle negative integers and zero themselves, if needed.
Math libraries tend to either handle it the same way, or restrict the input to positive numbers.
Wrapping up
I feel as though I had to walk a line, here, between too much and too little. I stopped short of doing all of the derivations to get to the formulas I used, or computing my own coefficients, in favor of a hopefully more digestible look at the topic. Please do check the External Resources, below, if you’d really like to get down in the weeds. And as always, please feel free to contact me with your feedback!
External Resources
- Gamma function at Wolfram Mathworld
- Lanczos approximation at Wikipedia
- The Lanczos Approximation at UCSF
- The Lanczos Approximation for the Γ-Function with Complex Coefficients by William Rea
- Complex Numbers at Math is Fun