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 _{k=0}^{\mathrm{\infty }}\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}^{\mathrm{\infty }}\frac{1}{{3888}^n}\frac{(1794n-297)}{n(2n-1)}\frac{n!(\frac{1}{2}{)}_n}{(\frac{1}{6}{)}_n(\frac{5}{6}{)}_n}.$$

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

$$\log (10)=46\mathrm{atanh}(1/31)+34\mathrm{atanh}(1/49)+20\mathrm{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}^{\mathrm{\infty }}\frac{H_k}{(k!{)}^2}{\left(\frac{x}{2}\right)}^{2k},I_0(x)=\sum _{k=0}^{\mathrm{\infty }}\frac{1}{(k!{)}^2}{\left(\frac{x}{2}\right)}^{2k}$$

$$2x{I}_0(x)K_0(x)\sim \sum _{k=0}^{\mathrm{\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}^{\mathrm{\infty }}\frac{(-4096{)}^{k}P(k)}{k^3(2k-1)(3k-1)(3k-2)(6k-1)(6k-5)(\genfrac{}{}{0ex}{}{5k}{k})(\genfrac{}{}{0ex}{}{10k}{5k})(\genfrac{}{}{0ex}{}{12k}{6k})}$$

where

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

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}^{\mathrm{\infty }}\frac{(-1{)}^{k-1}P(k)}{k^5(2k-1{)}^3(3k-1)(3k-2)(4k-1)(4k-3)(6k-1)(6k-5)(\genfrac{}{}{0ex}{}{5k}{k})(\genfrac{}{}{0ex}{}{5k}{2k})(\genfrac{}{}{0ex}{}{9k}{4k})(\genfrac{}{}{0ex}{}{10k}{5k})(\genfrac{}{}{0ex}{}{12k}{6k})}$$

where

$$\begin{array}{r}\begin{array}{cc}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{array}\end{array}$$

discovered by Zuniga [Zun2023].

Khinchin’s constant

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

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

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{array}{r}\begin{array}{c}\sum _{n=M+1}^{\mathrm{\infty }}\zeta (2n,N)=\sum _{n=M+1}^{\mathrm{\infty }}\sum _{k=0}^{\mathrm{\infty }}(k+N{)}^{-2n}\le \sum _{n=M+1}^{\mathrm{\infty }}(N^{-2n}+{\int }_0^{\mathrm{\infty }}(t+N{)}^{-2n}dt)\\ =\sum _{n=M+1}^{\mathrm{\infty }}\frac{1}{N^{2n}}(1+\frac{N}{2n-1})\le \sum _{n=M+1}^{\mathrm{\infty }}\frac{N+1}{N^{2n}}=\frac{1}{N^{2M}(N-1)}\le \frac{1}{N^{2M}}.\end{array}\end{array}$$

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 }^{\prime }(-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}^{\mathrm{\infty }}\frac{1}{F_n}=\sum _{n=0}^{\mathrm{\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/\varphi {)}^{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.