# 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

$\Gamma(1/3) = \left( \frac{12 \pi^4}{\sqrt{10}} \sum_{k=0}^{\infty} \frac{(6k)!(-1)^k}{(k!)^3 (3k)! 3^k 160^{3k}} \right)^{1/6}, \quad \quad \Gamma(1/4) = \sqrt{\frac{(2\pi)^{3/2}}{\operatorname{agm}(1, \sqrt 2)}}.$

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.