Algorithms for the gamma function

The Stirling series

In general, the gamma function is computed via the Stirling series

\[\log \Gamma(z) = \left(z-\frac{1}{2}\right)\log z - z + \frac{\ln {2 \pi}}{2} + \sum_{k=1}^{n-1} \frac{B_{2k}}{2k(2k-1)z^{2k-1}} + R(n,z)\]

where ([Olv1997] pp. 293-295) the remainder term is exactly

\[R_n(z) = \int_0^{\infty} \frac{B_{2n} - {\tilde B}_{2n}(x)}{2n(x+z)^{2n}} dx.\]

To evaluate the gamma function of a power series argument, we substitute \(z \to z + t \in \mathbb{C}[[t]]\).

Using the bound for \(|x+z|\) given by [Olv1997] and the fact that the numerator of the integrand is bounded in absolute value by \(2 |B_{2n}|\), the remainder can be shown to satisfy the bound

\[|[t^k] R_n(z+t)| \le 2 |B_{2n}| \frac{\Gamma(2n+k-1)}{\Gamma(k+1) \Gamma(2n+1)} \; |z| \; \left(\frac{b}{|z|}\right)^{2n+k}\]

where \(b = 1/\cos(\operatorname{arg}(z)/2)\). Note that by trigonometric identities, assuming that \(z = x+yi\), we have \(b = \sqrt{1 + u^2}\) where

\[u = \frac{y}{\sqrt{x^2 + y^2} + x} = \frac{\sqrt{x^2 + y^2} - x}{y}.\]

To use the Stirling series at \(p\)-bit precision, we select parameters \(r\), \(n\) such that the remainder \(R(n,z)\) approximately is bounded by \(2^{-p}\). If \(|z|\) is too small for the Stirling series to give sufficient accuracy directly, we first translate to \(z + r\) using the formula \(\Gamma(z) = \Gamma(z+r) / (z (z+1) (z+2) \cdots (z+r-1))\).

To obtain a remainder smaller than \(2^{-p}\), we must choose an \(r\) such that, in the real case, \(z + r > \beta p\), where \(\beta > \log(2) / (2 \pi) \approx 0.11\). In practice, a slightly larger factor \(\beta \approx 0.2\) more closely balances \(n\) and \(r\). A much larger \(\beta\) (e.g. \(\beta = 1\)) could be used to reduce the number of Bernoulli numbers that have to be precomputed, at the expense of slower repeated evaluation.

Rational arguments

We use efficient methods to compute \(y = \Gamma(p/q)\) where \(q\) is one of \(1, 2, 3, 4, 6\) and \(p\) is a small integer.

The cases \(\Gamma(1) = 1\) and \(\Gamma(1/2) = \sqrt \pi\) are trivial. We reduce all remaining cases to \(\Gamma(1/3)\) or \(\Gamma(1/4)\) using the following relations:

\[\Gamma(2/3) = \frac{2 \pi}{3^{1/2} \Gamma(1/3)}, \quad \quad \Gamma(3/4) = \frac{2^{1/2} \pi}{\Gamma(1/4)},\]

\[\Gamma(1/6) = \frac{\Gamma(1/3)^2}{(\pi/3)^{1/2} 2^{1/3}}, \quad \quad \Gamma(5/6) = \frac{2 \pi (\pi/3)^{1/2} 2^{1/3}}{\Gamma(1/3)^2}.\]

We compute \(\Gamma(1/3)\) and \(\Gamma(1/4)\) rapidly to high precision using

\[\frac{\pi}{\Gamma(1/3)^3} = \frac{1}{960 \cdot 10^{1/4}} \sum_{k=0}^{\infty} \left(-\frac{3^2}{2^9 \cdot 5^3}\right)^k (9108 k + 279) \frac{(1/12)_k (5/12)_k}{(k!)^2}, \quad \quad \Gamma(1/4) = \sqrt{\frac{(2\pi)^{3/2}}{\operatorname{agm}(1, \sqrt 2)}}.\]

where the infinite series is due to Guillera. An alternative formula which could be used for \(\Gamma(1/3)\) is

\[\Gamma(1/3) = \frac{2^{4/9} \pi^{2/3}}{3^{1/12} \left( \operatorname{agm}\left(1,\frac{1}{2} \sqrt{2+\sqrt{3}}\right)\right)^{1/3}},\]

but this appears to be slightly slower in practice.