.. _algorithms_constants:

Algorithms for mathematical constants
===============================================================================

Most mathematical constants are evaluated using the generic hypergeometric
summation code.

Pi
-------------------------------------------------------------------------------

`\pi` is computed using the Chudnovsky series

    .. math ::

        \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

.. math ::

    \log(2) = \frac{3}{4} \sum_{k=0}^{\infty} \frac{(-1)^k (k!)^2}{2^k (2k+1)!}

.. math ::

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

.. math ::

    \gamma = \frac{S_0(2n) - K_0(2n)}{I_0(2n)} - \log(n)

in which `n` is a free parameter and

.. math ::

    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}

.. math ::

    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

.. math ::

    C = \frac{1}{64} \sum_{k=1}^{\infty} \frac{256^k (580k^2-184k+15)}{k^3(2k-1){6k\choose 3k}{6k\choose 4k}{4k\choose 2k}}

given in [PP2010]_.

Khinchin's constant
-------------------------------------------------------------------------------

Khinchin's constant `K_0` is computed using the formula

.. math ::

    \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

.. math ::

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

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

.. math ::

    \zeta(3) = \frac{1}{64} \sum_{k=0}^\infty
        (-1)^k (205k^2 + 250k + 77) \frac{(k!)^{10}}{[(2k+1)!]^5}.