# 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 development1 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 would be either taking the logarithmic derivative of either of the Gamma function approximations, or 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.

1. Yang ZH, Tian JF. An accurate approximation formula for gamma function. J Inequal Appl. 2018;2018(1):56. doi: 10.1186/s13660-018-1646-6. Epub 2018 Mar 6. PMID: 29540975; PMCID: PMC5840229. ↩︎