Mathematical constants

Most of the mathematical constants provided by arb are sums of hypergeometric series, series whose term ratio is a rational function of the index, and are evaluated by binary splitting. The series for $\pi$, $e$, $\log 2$, $\log 10$, $\gamma$, Catalan's constant, $\zeta(3)$ and $\zeta(5)$ are summed by the mp_real module, and the arb_const_ functions call it through _mp_real_const_arb; the formulas were chosen as the fastest among those collected by y-cruncher [Yee2025]. The implementation side (splitting trees, threads, the shared cache of exact floors) is described in Constants in mp_real; the sections below give the mathematics and the remaining constants.

Binary splitting#

Consider a sum $S = \sum_{k=0}^{N-1} a(k)/b(k) \cdot \prod_{j<k} p(j)/q(j)$ with polynomials $a, b, p, q$ with integer coefficients. For $0 \le u < v \le N$ define $P(u, v) = \prod_{j=u}^{v-1} p(j)$, $Q(u, v) = \prod_{j=u}^{v-1} q(j)$, and $T(u, v) = Q(u, v) \sum_{k=u}^{v-1} (\ldots)$, the partial sum scaled to an integer. Then for $u < w < v$

(with the obvious extension for $a/b$), so the exact sum is obtained by a balanced recursion and one final division $T/Q$.

The series below add between 2 and 47 bits per term. The constant factor is what distinguishes the formulas, and it is governed by the number of bits per term and the size of the coefficients [HP1998] [BZ2011].

Caching#

Each arb_const_ function keeps a thread-local cached ball. A request at precision $p$ above the cached precision recomputes the constant at $p + 32$ bits; the result is then rounded to $p$ bits, and later requests at precision at most $p$ only round the cached value. $\pi$, $\log 2$ and $\gamma$ are moreover read from static tables (of $\pi/4$ and $\log 2$ to 4608 bits, of $\gamma$ to 3456 bits) below these sizes, correctly rounded to nearest with a radius of half a unit in the last place. The caches are freed by flint_cleanup.

$\pi$#

$\pi$ is computed by the Chudnovsky series [CC1988]

which adds about 47 bits (14 digits) per term. In the hypergeometric form used by mp_real, $\pi = S / \sqrt{10005}$ with $R(k) = (6k - 5)(2k - 1)(6k - 1)$, $P(k) = (13591409 + 545140134 k) R(k)$ and $Q(k) = -(640320^3/24) k^3$. Methods based on the arithmetic-geometric mean [Sal1976] [Bre1976agm] have a better asymptotic complexity, $O(M(n) \log n)$, but are about three times slower in practice. arb_const_sqrt_pi and arb_const_log_sqrt2pi are computed from arb_const_pi with arb_sqrt and arb_log.

$e$#

$e = \sum_{k \ge 0} 1/k!$ is summed by binary splitting, with $P = 1$ and $Q(k) = k$; this is the simplest hypergeometric series and its cost is dominated by the final division.

Logarithms of integers#

$\log 2$ is computed by Zuniga's series [Zun2025]

which adds 11.9 bits per term. (Earlier versions used another series of Zuniga [Zun2023b], with terms decreasing like $3888^{-k}$.) $\log 10 = \log 2 + \log 5$ uses the Machin-type formulas for the first three primes in terms of

for example $\log 10 = 46 \operatorname{atanh}(1/31) + 34 \operatorname{atanh}(1/49) + 20 \operatorname{atanh}(1/161)$, with the three series evaluated in parallel. The logarithms of more primes, used by the elementary functions (Elementary functions), come from larger sets of such formulas (arb_log_primes_vec_bsplit).

Euler's constant#

Euler's constant $\gamma$ is computed by the Brent-McMillan formula [BM1980] [MPFR2012]

where $n$ is a free parameter and

The first two series are summed together up to $k = N - 1$ with $N \ge \alpha n + 1$, where $\alpha \approx 4.9706257595442318644$ satisfies $\alpha(\log \alpha - 1) = 3$, and the asymptotic series up to $k = 2n - 1$. With these parameters the error is bounded by $24 e^{-8n}$ [BJ2013]. The current implementation sums $S_0$ and $I_0$ in one binary splitting over dual numbers, and chooses $n$ smooth over $\{2, 3, 5, 7\}$ so that $\log n$ is cheap; see Constants in mp_real.

Catalan's constant#

Catalan's constant is computed by the series of Pilehrood and Pilehrood [PP2010]

which adds about 5.4 bits per term and was measured to be the fastest available formula.

$\zeta(3)$ and $\zeta(5)$#

Apéry's constant $\zeta(3)$ is computed by a hypergeometric series found by Zuniga [Zun2023],

with $P(k)$ a polynomial of degree 11 given in the source. $\zeta(5)$ is computed from Z.-W. Sun's identity

with the three parts $S$, $\zeta(3)$ and $\pi$ computed in parallel when threads are available.

Khinchin's constant#

Khinchin's constant $K_0$ is computed from the formula [BBC1997]

where $N \ge 2$ is a free parameter (arb_const_khinchin uses $N \approx p^{0.35}$). The finite sum has a minus sign: each product of logarithms is negative.

For an error of $2^{-p}$ it therefore suffices to take $M \ge p / (2 \log_2 N)$. The Hurwitz zeta values are computed as $\zeta(2n, N) = \zeta(2n) - \sum_{k=1}^{N-1} k^{-2n}$, with arb_zeta_ui and a vector of powers $k^{-2n}$ that is updated from one $n$ to the next.

Glaisher's constant#

The Glaisher-Kinkelin constant $A = \exp(1/12 - \zeta'(-1))$ is computed directly from this formula, with $\zeta'(-1)$ from acb_dirichlet_zeta_jet (see Zeta functions, polylogarithms and L-functions). Up to a cutoff (_acb_dirichlet_l_jet_use_afe: $3.5 \cdot 10^5$ bits with one thread, about $8 \cdot 10^4$ bits with 8 threads), this uses Euler-Maclaurin summation directly at $s = -1$; the reflection formula is not used for $-4 < \operatorname{Re}(s) < 0$, since the arithmetic of Euler-Maclaurin summation is cheaper at $s = -1$ than at $s = 2$ (by a factor 1.3–1.7 for $\zeta'(-1)$).

Euler-Maclaurin summation needs $O(p)$ terms and $O(p)$ Bernoulli numbers at precision $p$, which costs $\tilde O(p^2)$ time and $O(p^2)$ bits of memory (several gigabytes at $p = 4 \cdot 10^5$). Above the cutoff, acb_dirichlet_zeta_jet instead computes $\zeta(-1+x)$ to $O(x^2)$ with the approximate functional equation (acb_dirichlet_l_fmpq_afe_jet), a sum of $O(\sqrt{p})$ incomplete gamma function jets evaluated by binary splitting, in time $\tilde O(p^{1.5})$ and with $O(p)$ memory.

The reciprocal Fibonacci constant#

$\sum_{n \ge 1} 1/F_n$ is computed from Gosper's series [Gos1974], corrected in [Arn2012],

where $L_n = 2F_{n-1} + F_n$ is a Lucas number. The truncation error after $N \ge 1$ terms is bounded by $(1/\phi)^{N^2}$, so $N \approx 1.2 \sqrt{p}$ terms give $p$ bits. The series is not hypergeometric, but factoring out the multiplicative recurrence for $L_1 L_3 \cdots$ turns it into a product of $O(\sqrt{p})$ matrices with $O(\sqrt{p})$-bit entries, which binary splitting evaluates in quasilinear time.

Use in FLINT#

The complex versions in acb (for example acb_const_pi) call these functions. The constants $\Gamma(1/3)$ and $\Gamma(1/4)$, used by arb_hypgeom_gamma_fmpq, are also computed by mp_real; see Constants in mp_real and the gamma function at rational arguments.

The tutorial Computing π and other constants computes $\pi$ and other constants to many digits in complete programs.

@misc{flint-doc-algorithms-special-functions-constants,
  author       = {{The FLINT team}},
  title        = {{Mathematical constants}},
  howpublished = {FLINT: Fast Library for Number Theory, documentation, version 3.7.0-dev},
  year         = {2026},
  url          = {https://flintlib.org/doc/algorithms/special-functions/constants.html},
  note         = {Generated 2026-11-02}
}