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}.\]