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 + t^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.
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
An alternative formula which could be used for \(\Gamma(1/3)\) is
but this appears to be slightly slower in practice.