Mathematical constants

mp_real computes twelve constants: $\pi/4$, $\log 2$, Euler's constant $\gamma$, $e$, $\log 10$, Catalan's constant $G$, $\zeta(3)$, $\zeta(5)$, $\Gamma(1/3)$, $\Gamma(1/4)$, $2/\pi$ and $1/\log 2 - 1$. They are the implementations behind arb_const_pi, arb_const_log2, arb_const_euler, arb_const_e, arb_const_log10, arb_const_catalan, arb_const_apery, arb_const_zeta5 and the $\Gamma$ values of arb_hypgeom_gamma_fmpq, and they are used inside the module by the argument reductions of the elementary functions.

Two forms#

The ball functions mp_real_const_* return a ball for the constant accurate to about $n$ limbs. The limb functions _mp_real_const_* return the exact floor $\lfloor c B^n \rfloor$ as a fixed-point number with $n$ fraction limbs, plus a units limb for the constants above 1, and set err to 1: the true value lies in $[y, y + 1)$ ulps. An exact floor is a canonical value, independent of how it was computed. It is determined from a ball computed with a few guard limbs, recomputed with more guard limbs in the rare case that the ball contains a multiple of $B^{-n}$.

Caching#

With cache = 0 the constant is computed from scratch and nothing is kept. With cache = 1 the value comes from a process-wide cache that holds the floor at the largest precision requested so far. Floors nest: the top limbs of $\lfloor cB^N \rfloor$ are $\lfloor cB^n \rfloor$ for $n \le N$, so a shorter request is a copy of the top limbs. A longer request recomputes the entry at $\max(n + 5, 1.5N)$ limbs, so that a sequence of growing requests costs a bounded factor more than the largest. All threads read the same entry; one that needs a longer entry builds it under the cache's lock, and readers never wait. The ball version with cache = 1 returns the cached floor with a radius of one ulp. $\pi/4$, $\log 2$, $2/\pi$ and $1/\log 2 - 1$ also come from static tables of 64 limbs, so that the elementary functions at small precision touch no cache at all. _mp_real_const_clear_cache (or flint_cleanup_master) frees the caches.

The arb constants call _mp_real_const_arb, which computes with cache = 0 and leaves caching to arb's own caches.

Algorithms#

All constants except $\log 10$ are hypergeometric series summed by mp_real_hypgeom_series, with formulas chosen from y-cruncher's formula collection as the fastest in this implementation:

The independent parts (the series and $\pi$ in $\zeta(5)$ and the $\Gamma$ values, the three logarithms of $\log 10$, the logarithms needed by $\gamma$) run on separate threads. The details, including the dual-number splitting for $\gamma$, are in Constants in mp_real.

Example#

#include <flint/flint.h>
#include <flint/mp_real.h>
#include <flint/arb.h>

int main(void)
{
    mp_real_t c;
    arb_t a;
    ulong f[3], err;
    slong n = mp_real_prec_bits(200);

    mp_real_init(c);
    arb_init(a);

    /* zeta(3) as a ball, computed from scratch */
    mp_real_const_zeta3(c, n, 0);
    mp_real_get_arb(a, c);
    flint_printf("zeta(3) = ");
    arb_printn(a, 50, 0);
    flint_printf("\n");

    /* the cached floor of e B^2: two fraction limbs and a units limb */
    _mp_real_const_e(f, &err, 2, 1);
    flint_printf("floor(e B^2) = %wu %016wx %016wx, err = %wu\n",
        f[2], f[1], f[0], err);

    /* Catalan's constant through the arb entry point */
    _mp_real_const_arb(a, mp_real_const_catalan, 128);
    flint_printf("G = ");
    arb_printn(a, 30, 0);
    flint_printf("\n");

    mp_real_clear(c);
    arb_clear(a);
    flint_cleanup_master();
    return 0;
}

zeta(3) = [1.2020569031595942853997381615114499907649862923405 +/- 1.12e-51]
floor(e B^2) = 2 b7e151628aed2a6a bf7158809cf4f3c7, err = 1
G = [0.915965594177219015054603514932 +/- 3.85e-31]

The example assumes 64-bit limbs; the hexadecimal digits of $e$ are $2.\mathrm{B7E151628AED2A6ABF7158809CF4F3C7}\ldots$