# 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¶

The standalone constant $$\log(2)$$ is computed using Zuniga’s series [Zun2023b]

$\log(2) = \frac{1}{2} \sum_{n=1}^\infty \frac{1}{3888^n} \frac{(1794 n-297)}{n(2n-1)} \frac{n! (\tfrac{1}{2})_n}{(\tfrac{1}{6})_n (\tfrac{5}{6})_n}.$

Logarithms of other small integers are in certain situations computed using Machin-like formulas, e.g.:

$\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 = \frac{1}{768} \sum_{k=1}^{\infty} \frac{(-4096)^k P(k)} {k^3 (2k-1)(3k-1)(3k-2)(6k-1)(6k-5) {5k \choose k} {10k \choose 5k} {12k \choose 6k}}$

where

$\begin{split}\begin{matrix} P(k) & = -43203456k^6 + 92809152k^5 - 76613904k^4 \\ & + 30494304k^3 - 6004944k^2 + 536620^k - 17325, \end{matrix}\end{split}$

discovered by Zuniga [Zun2023]. It was previously computed using a series given in [PP2010].

## Apery’s constant¶

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

$\zeta(3) = \frac{1}{48} \sum_{k=1}^{\infty} \frac{(-1)^{k-1} P(k)}{k^5 (2k-1)^3(3k-1)(3k-2)(4k-1)(4k-3)(6k-1)(6k-5){5k \choose k}{5k \choose 2k}{9k \choose 4k}{10k \choose 5k}{12k \choose 6k}}$

where

$\begin{split}\begin{matrix} P(k) & = 1565994397644288k^{11} - 6719460725627136k^{10} + 12632254526031264k^9 \\ & - 13684352515879536k^8 + 9451223531851808k^7 - 4348596587040104k^6 \\ & + 1352700034136826k^5 - 282805786014979k^4 + 38721705264979k^3 \\ & - 3292502315430k^2 + 156286859400k - 3143448000, \end{matrix}\end{split}$

discovered by Zuniga [Zun2023].

## 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$$.

## Reciprocal Fibonacci constant¶

We use Gosper’s series ([Gos1974], corrected in [Arn2012])

$\sum_{n=1}^{\infty} \frac{1}{F_n} = \sum_{n=0}^{\infty} \frac{(-1)^{n(n-1)/2} (F_{4n+3} + (-1)^n F_{2n+2})}{F_{2n+1} F_{2n+2} L_1 L_3 \cdots L_{2n+1}}$

where $$L_n = 2F_{n-1} + F_n$$ denotes a Lucas number. The truncation error after $$N \ge 1$$ terms is bounded by $$(1 / \phi)^{N^2}$$. The series is not of hypergeometric type, but we can evaluate it in quasilinar time using binary splitting; factoring out a multiplicative recurrence for $$L_1 L_3 \cdots$$ allows computing the series as a product of $$O(\sqrt{p})$$ matrices with $$O(\sqrt{p})$$-bit entries.