# Algorithms for mathematical constants¶

Most mathematical constants are evaluated using the generic hypergeometric summation code.

## Pi¶

$$\pi$$ is computed using the Chudnovsky series

$\frac{1}{\pi} = 12 \sum^\infty_{k=0} \frac{(-1)^k (6k)! (13591409 + 545140134k)}{(3k)!(k!)^3 640320^{3k + 3/2}}$

which is hypergeometric and adds roughly 14 digits per term. Methods based on the arithmetic-geometric mean seem to be slower by a factor three in practice.

A small trick is to compute $$1/\sqrt{640320}$$ instead of $$\sqrt{640320}$$ at the end.

## Logarithms of integers¶

We use the formulas

$\log(2) = \frac{3}{4} \sum_{k=0}^{\infty} \frac{(-1)^k (k!)^2}{2^k (2k+1)!}$
$\log(10) = 46 \operatorname{atanh}(1/31) + 34 \operatorname{atanh}(1/49) + 20 \operatorname{atanh}(1/161)$

## Euler’s constant¶

Euler’s constant $$\gamma$$ is computed using the Brent-McMillan formula ([BM1980], [MPFR2012])

$\gamma = \frac{S_0(2n) - K_0(2n)}{I_0(2n)} - \log(n)$

in which $$n$$ is a free parameter and

$S_0(x) = \sum_{k=0}^{\infty} \frac{H_k}{(k!)^2} \left(\frac{x}{2}\right)^{2k}, \quad I_0(x) = \sum_{k=0}^{\infty} \frac{1}{(k!)^2} \left(\frac{x}{2}\right)^{2k}$
$2x I_0(x) K_0(x) \sim \sum_{k=0}^{\infty} \frac{[(2k)!]^3}{(k!)^4 8^{2k} x^{2k}}.$

All series are evaluated using binary splitting. The first two series are evaluated simultaneously, with the summation taken up to $$k = N - 1$$ inclusive where $$N \ge \alpha n + 1$$ and $$\alpha \approx 4.9706257595442318644$$ satisfies $$\alpha (\log \alpha - 1) = 3$$. The third series is taken up to $$k = 2n-1$$ inclusive. With these parameters, it is shown in [BJ2013] that the error is bounded by $$24e^{-8n}$$.

## Catalan’s constant¶

Catalan’s constant is computed using the hypergeometric series

$C = \sum_{k=0}^{\infty} \frac{(-1)^k 4^{4 k+1} \left(40 k^2+56 k+19\right) [(k+1)!]^2 [(2k+2)!]^3}{(k+1)^3 (2 k+1) [(4k+4)!]^2}$

## Khinchin’s constant¶

Khinchin’s constant $$K_0$$ is computed using the formula

$\log K_0 = \frac{1}{\log 2} \left[ \sum_{k=2}^{N-1} \log \left(\frac{k-1}{k} \right) \log \left(\frac{k+1}{k} \right) + \sum_{n=1}^\infty \frac {\zeta (2n,N)}{n} \sum_{k=1}^{2n-1} \frac{(-1)^{k+1}}{k} \right]$

where $$N \ge 2$$ is a free parameter that can be used for tuning [BBC1997]. If the infinite series is truncated after $$n = M$$, the remainder is smaller in absolute value than

\begin{align}\begin{aligned}\sum_{n=M+1}^{\infty} \zeta(2n, N) = \sum_{n=M+1}^{\infty} \sum_{k=0}^{\infty} (k+N)^{-2n} \le \sum_{n=M+1}^{\infty} \left( N^{-2n} + \int_0^{\infty} (t+N)^{-2n} dt \right)\\= \sum_{n=M+1}^{\infty} \frac{1}{N^{2n}} \left(1 + \frac{N}{2n-1}\right) \le \sum_{n=M+1}^{\infty} \frac{N+1}{N^{2n}} = \frac{1}{N^{2M} (N-1)} \le \frac{1}{N^{2M}}.\end{aligned}\end{align}

Thus, for an error of at most $$2^{-p}$$ in the series, it is sufficient to choose $$M \ge p / (2 \log_2 N)$$.

## Glaisher’s constant¶

Glaisher’s constant $$A = \exp(1/12 - \zeta'(-1))$$ is computed directly from this formula. We don’t use the reflection formula for the zeta function, as the arithmetic in Euler-Maclaurin summation is faster at $$s = -1$$ than at $$s = 2$$.

## Apery’s constant¶

Apery’s constant $$\zeta(3)$$ is computed using the hypergeometric series

$\zeta(3) = \frac{1}{64} \sum_{k=0}^\infty (-1)^k (205k^2 + 250k + 77) \frac{(k!)^{10}}{[(2k+1)!]^5}.$