Introduction⌗

The factorial function and its extension to the complex plane - the Gamma function are without doubt my favourite mathematical devices. In this blog post, I would like to shed some light on approximation methods for the Gamma function, as well as functions closely related to it (polygamma, digamma, etc…).

The Gamma function⌗

The Gamma function is canonically defined using an integral as follows:

$$\Gamma(z) = \int_0^\infty t^{z-1}e^{-t}dt$$

The Gamma function is a generalization of the factorial function, therefore the following holds:

$$\Gamma(n+1) = n!$$

Since $\Gamma(1) = 1$ and $\Gamma(n+1)=n \Gamma(n)$, which follows from the integral definition.

Euler and Weierstrass respectively suggest the following infinite product expansions for the Gamma function (where $\gamma$ is the Euler-Mascheroni constant), however they’re not that useful in practice:

$${\displaystyle \Gamma (z)={\frac {1}{z}}\prod _{n=1}^{\infty }{\frac {\left(1+{\frac {1}{n}}\right)^{z}}{1+{\frac {z}{n}}}} = {\frac {e^{-\gamma z}}{z}}\prod _{n=1}^{\infty }\left(1+{\frac {z}{n}}\right)^{-1}e^{z/n}}$$

A more useful tool is the reflection formula:

$${\displaystyle \Gamma (1-z)\Gamma (z)={\frac {\pi }{\sin \pi z}}}$$

And consequently:

$$\Gamma (-z) = -\frac{\Gamma(1 - z)}{z} = -\frac {\pi }{z \Gamma(z) \sin \pi z} = - \frac{\pi \csc(\pi z)}{\Gamma(z + 1)}$$

Which lets us assume that $Re(z) > 0$ without loss of generality.

The Digamma function⌗

The digamma function is a logarithmic derivative of the Gamma function, and is consequently defined as follows:

$$\psi(z) = \frac{d}{dz} \ln \Gamma(z) = \frac{\Gamma’(z)}{\Gamma(z)}$$

The digamma function has its own reflection and recurrence formulae:

$${\displaystyle \psi (1-x)-\psi (x)=\pi \cot \pi x}$$ $${\displaystyle \psi (x+1)=\psi (x)+{\frac {1}{x}}}$$

The Polygamma function⌗

The $m$-th polygamma function is defined as the $m$-th derivative of the digamma function:

$${\displaystyle \psi ^{(m)}(z)={\frac {\mathrm {d} ^{m}}{\mathrm {d} z^{m}}}\psi (z)}$$

An integral representation for the polygamma function assuming $m > 0$ and $Re(z) > 0$ is given by:

$${\displaystyle {\psi ^{(m)}(z)=(-1)^{m+1}\int _{0}^{\infty }{\frac {t^{m}e^{-zt}}{1-e^{-t}}}dt=-\int _{0}^{1}{\frac {t^{z-1}}{1-t}}(\ln t)^{m}dt}}$$

The recurrence relation and reflection formula follow:

$$\psi ^{{(m)}}(z+1)=\psi ^{{(m)}}(z)+{\frac {(-1)^{m}\ m!}{z^{{m+1}}}}$$ $${\displaystyle (-1)^{m}\psi ^{(m)}(1-z)-\psi ^{(m)}(z)=\pi {\frac {\mathrm {d} ^{m}}{\mathrm {d} z^{m}}}\cot {\pi z}}$$

Stirling’s approximation⌗

Likely the most widely known formula for approximating the factorial function:

$${\displaystyle n!\sim {\sqrt {2\pi n}}\left({\frac {n}{e}}\right)^{n}}$$

Stirling approximation is only ever useful for large values of $n$. It can be used to approximate the gamma function for positive real argument, and as a consequence of the reflection formula, negative real argument:

$$\Gamma(-z) = - \frac{\pi \csc(\pi z)}{{\sqrt {2\pi (z + 1)}}\left({\frac {z+1}{e}}\right)^{z+1}}$$

Because of how inaccurate the approximation is for smaller values, one can use the recurrence relation obtain a large argument.

Ramanjuan’s formula⌗

A less commonly known formula derived by Ramanujan as a refinement to the Stirling’s approximation:

$$n! \sim \sqrt{\pi} \left(\frac{n}{e}\right)^n \sqrt[6]{8n^3 + 4n^2 + n + \frac{1}{30}}$$

Ramanujan’s formula is more accurate than the Stirling’s formula, especially for smaller values of $n$ (e.g. $n = 5$).

Yang & Tian’s formula⌗

A very recent development¹ by Zhen-Hang Yang and Jing-Feng Tian, argued to be even more accurate than Ramanujan’s formula.

$$\Gamma\left(x+1\right)\sim\sqrt{2 \pi x}\left(\frac{x}{e}\right)^{x}\left(x \sinh \frac{1}{x}\right)^{\frac{x}{2}}\exp\left(\frac{7}{324 x^{3} \left( 35 x^{2} + 33 \right)}\right)$$

Spouge approximation⌗

Spouge’s approximation is a modification of Stirling’s approximation. It’s, in my opinion, the most efficient way to compute the values of the Gamma function when using arbitrary precision arithmetics. The formula follows:

$$\Gamma (z+1)=(z+a)^{z+{\frac {1}{2}}}e^{-z-a}\left(c_{0}+\sum_{k=1}^{a-1}{\frac {c_{k}}{z+k}}\right)$$

The coefficients are defined as:

$${\displaystyle {\begin{aligned}c_{0}&={\sqrt {2\pi }}\\c_{k}&={\frac {(-1)^{k-1}}{(k-1)!}}(-k+a)^{k-{\frac {1}{2}}}e^{-k+a}\qquad k\in {1,2,\dots ,a-1}.\end{aligned}}}$$

Approximating the gamma function for small arguments⌗

The following modification of Stirling’s formula gives satisfactory results near $x = 0$:

$$\Gamma(x)\approx\sqrt{2\pi\left(x+\tfrac16\right)}\ \frac{x^{x-1}}{e^x}$$

Approximating the digamma function⌗

My preferred method of approximating the digamma function for small arguments would be using McCullagh’s series. As an example of the first method, we will approximate $\Gamma(x)$ using the Ramanujan approximation for $(x - 1)!$

$$\begin{aligned} \psi(n) &\sim \frac{d}{dn} \ln \sqrt{\pi} \left(\frac{n-1}{e}\right)^{n-1} \sqrt[6]{8 n^3 - 20n^2 + 17n - \frac{149}{30}}\\ &= \frac{40 n (3 n - 5) + 85}{30 n (8n^2 - 20n + 17) - 149} + \log(n-1) \end{aligned}$$

There are also simpler to compute asymptotic approximations for the digamma function, however, they’re not particularly useful because of their inaccuracy. For example:

$$\psi(x) \sim \ln x - \frac{1}{2x}$$

Finally, the logarithmic derivative of the Spouge’s approximation also performs well. Temporarily fix $f(z)$ as the coefficient term:

$$\begin{aligned}\psi(z+1) &\sim \frac{d}{dz} \ln((z+a)^{z + \frac{1}{2}}e^{-z-a}f(z))\\ &= \frac{1-2a}{2a+2z} + \log(a+z) + \frac{f’(z)}{f(z)}\\ &= \frac{1-2a}{2a+2z} + \log(a+z) - \frac{\sum\limits_{k=1}^{a-1}{\frac {c_{k}}{(z+k)^2}}}{c_{0}+\sum\limits_{k=1}^{a-1}{\frac {c_{k}}{z+k}}} \end{aligned}$$

However, McCullagh’s series for the Digamma function let us solve this problem more elegantly:

$$\psi(z+1) = -\gamma - \sum_{n=1}^{\infty} (-z)^n \left(d_n + \frac{c_n}{z + n}\right)$$

Where the coefficients can be precomputed:

$$c_n = \frac{1}{n^n}\text{, }d_n=\zeta(n+1, n+1)=\sum_{k=n+1}^{\infty} \frac{1}{k^{n+1}}$$

The accuracy of McCullagh’s series is the best for $\psi(1)$ and dimnishes as the argument grows.

For larger arguments, I use the digamma reflection formula to increase the argument and then use the asymptotic expansion:

public static double digamma(double x) {
    double digamma = 0;

    // Reflection:
    if (x < 0) {
        digamma -= Math.PI / Math.tan(Math.PI * x);
        x = 1 - x;
    }

    // Small enough arguments:
    if (x > 0 && x <= 1e-5)
        return digamma - 0.577215664901532860606512090082 - 1 / x;

    // Increase argument:
    while (x < 49) {
        digamma -= 1 / x;
        x += 1;
    }

    // Use asymptotics
    final double inv = 1 / (x * x);
    digamma += Math.log(x) - 0.5 / x + inv * (-1d / 12 + inv * (1d / 120 - (1d / 252) * inv));

    return digamma;
}