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

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

## Euler’s constant¶

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

in which \(n\) is a free parameter and

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

where

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

where

discovered by Zuniga [Zun2023].

## Khinchin’s constant¶

Khinchin’s constant \(K_0\) is computed using the formula

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

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])

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.