# Algorithms for the gamma function¶

## The Stirling series¶

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

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

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

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

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:

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

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

but this appears to be slightly slower in practice.