# 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

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

## 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\).

## Apery’s constant¶

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