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:

Γ(z)=0tz1etdt\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:

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

Since Γ(1)=1\Gamma(1) = 1 and Γ(n+1)=nΓ(n)\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:

Γ(z)=1zn=1(1+1n)z1+zn=eγzzn=1(1+zn)1ez/n {\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:

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

And consequently:

Γ(z)=Γ(1z)z=πzΓ(z)sinπz=πcsc(πz)Γ(z+1) \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)>0Re(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:

ψ(z)=ddzlnΓ(z)=Γ(z)Γ(z) \psi(z) = \frac{d}{dz} \ln \Gamma(z) = \frac{\Gamma’(z)}{\Gamma(z)}

The digamma function has its own reflection and recurrence formulae:

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

The Polygamma function

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

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

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

ψ(m)(z)=(1)m+10tmezt1etdt=01tz11t(lnt)mdt {\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:

ψ(m)(z+1)=ψ(m)(z)+(1)m m!zm+1 \psi ^{{(m)}}(z+1)=\psi ^{{(m)}}(z)+{\frac {(-1)^{m}\ m!}{z^{{m+1}}}} (1)mψ(m)(1z)ψ(m)(z)=πdmdzmcotπz {\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:

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

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

Γ(z)=πcsc(πz)2π(z+1)(z+1e)z+1 \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!π(ne)n8n3+4n2+n+1306 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 nn (e.g. n=5n = 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.

Γ(x+1)2πx(xe)x(xsinh1x)x2exp(7324x3(35x2+33)) \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:

Γ(z+1)=(z+a)z+12eza(c0+k=1a1ckz+k)\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:

c0=2πck=(1)k1(k1)!(k+a)k12ek+ak1,2,,a1. {\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=0x = 0:

Γ(x)2π(x+16) xx1ex \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 Γ(x)\Gamma(x) using the Ramanujan approximation for (x1)!(x - 1)!

ψ(n)ddnlnπ(n1e)n18n320n2+17n149306=40n(3n5)+8530n(8n220n+17)149+log(n1)\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:

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

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

ψ(z+1)ddzln((z+a)z+12ezaf(z))=12a2a+2z+log(a+z)+f(z)f(z)=12a2a+2z+log(a+z)k=1a1ck(z+k)2c0+k=1a1ckz+k\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:

ψ(z+1)=γn=1(z)n(dn+cnz+n)\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:

cn=1nndn=ζ(n+1,n+1)=k=n+11kn+1c_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 ψ(1)\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;

  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. ↩︎