Bernoulli numbers

The Bernoulli numbers B_n are the rational numbers defined by x/(e^x - 1) = \sum_{n \ge 0} B_n x^n/n!. They vanish for odd n \ge 3, alternate in sign for even n, and grow quickly:

The numerator of B_n thus has about n \log_2 n bits; that of B_{10^6} has about 15.8 million. For the classical theory see [IR1990, Chapter 15] and [Coh2007b, Chapter 9].

The denominator#

The denominator is therefore found from the divisors of n (arith_bernoulli_number_denom tests the primes up to n + 1, with a table for n \le 178). It is small: at most \prod_{d \mid n} (d + 1), so it has at most \sigma_0(n) \log_2(n+1) bits, and usually far fewer. All methods below therefore compute the numerator N_n = d_n B_n, an integer.

The zeta function method#

Since d_n is known, it suffices to compute B_n numerically with an absolute error below 1/(2 d_n), that is, to about \log_2 |B_n| + \log_2 d_n + 1 significant bits. With p bits of precision, the formula above needs n! and \pi^n (cheap: binary splitting and a power) and \zeta(n), which is 1 + O(2^{-n}) and for which very few terms are needed:

so J \approx 2^{p/(n-1)} terms suffice, which is O(n/\log n) for p \approx n \log_2 n. Equivalently one evaluates the Euler product \zeta(n)^{-1} = \prod_{q \le J} (1 - q^{-n}) over primes, with O(n/\log^2 n) factors. FLINT's arb_bernoulli_ui_zeta uses the Euler product when n is small compared with the precision and an asymptotic expansion of \zeta(n) otherwise, in ball arithmetic. _bernoulli_fmpq_ui_zeta multiplies the ball by d_n and returns the unique integer it contains (increasing the precision in the unlikely event that the ball contains two). Each power q^{-n} costs O(M(n \log n) \log n), so the total is O(n M(n \log n)), about n^2 \log^{2+o(1)} n bit operations.

The reverse iterator#

To compute all of B_0, B_2, \ldots, B_n, evaluating each from scratch would cost n^3. Bloemen's observation [Blo2009] is that the power table can be shared: if the powers j^{-n} are stored in fixed point with p fractional bits, then j^{-(n-2)} = j^{-n} \cdot j^2 is one multiplication by a word. Going down in n thus updates the table with cheap exact multiplications (going up would need divisions), and the precision and the number of terms needed decrease along the way.

Each step costs O(J) multiplications of a p-bit number by a word, O(n/\log n \cdot n \log n) = O(n^2) bit operations, so all the even Bernoulli numbers up to n cost O(n^3) in the worst case but with a very small constant, and blocks of a range are independent, which makes the computation easy to parallelise. bernoulli_fmpq_vec_no_cache splits the range into blocks, runs one reverse iterator per block (largest blocks first) on the threads of flint_set_num_threads, and the cache used by the arb special functions is filled this way.

Harvey's multimodular algorithm#

Harvey [Har2010] computes B_n \bmod p for a prime p in time O(p \log^{1+o(1)} p), independently of n:

1. If (p - 1) \mid n, p divides the denominator. Otherwise, by Kummer's congruence B_n/n \equiv B_m/m \pmod p for m = n \bmod (p-1), so it suffices to compute B_m/m \bmod p with 2 \le m \le p - 3.
2. Voronoi's congruence expresses (c^m - 1) B_m / m, for any c coprime to p, as a sum over x \in (\mathbb{Z}/p\mathbb{Z})^\times of x^{m-1} times a simple function of \lfloor cx/p \rfloor. Enumerating x as the powers of a primitive root g, and using the symmetry x \leftrightarrow p - x, gives a sum of (p-1)/2 terms, each one multiplication modulo p.
3. With c = 2 (when 2^m \not\equiv 1 \pmod p), the floor terms are the bits of the binary expansion of x/p. Harvey reorganises the sum so that the inner loop reads the expansion of 1/p word by word and accumulates into tables indexed by 8-bit chunks; the tables are combined with precomputed weights at the end. This makes the constant very small. If 2^m \equiv 1, the slower generic sum with c = g is used.

With all primes 5 \le p \le X, (p-1) \nmid n, up to X \approx n \log n, the product exceeds |N_n|, and the Chinese remainder theorem gives the numerator: O(n^2 \log^{O(1)} n) operations in total, but fully parallel and with smaller constants than the zeta method for large n. FLINT combines the two methods: the residues modulo a product M of primes give the low part r = N_n \bmod M, and the zeta method only needs to find the integer t = (N_n - r)/M, a number with a fraction 1 - \alpha of the bits of N_n. The default \alpha grows from 0 at n = 18000 to 0.28 for large n.

Asymptotically faster methods exist: Harvey's later algorithm computes B_n in O(n^{4/3 + o(1)}) bit operations [HarveyBern2014], and Brent and Harvey compute all of B_0, \ldots, B_n in quasi-optimal time with power series [BrentHarvey2013]; they are not implemented.

Short vectors#

For a vector B_0, \ldots, B_{n-1} with moderate n, arith uses:

• Ramanujan's recursion (_arith_bernoulli_number_vec_recursive, n < 700). Ramanujan's identity [Ram1911]

  \binom{m+3}{m} B_m = \begin{cases} \frac{m+3}{3} \\ -\frac{m+3}{6} \end{cases} - \sum_{j=1}^{\lfloor m/6 \rfloor} \binom{m+3}{m - 6j} B_{m-6j}, \quad m \equiv 0, 2 \text{ or } 4 \pmod 6,

  (with (m+3)/3 for m \equiv 0, 2 and -(m+3)/6 for m \equiv 4) needs only a sixth of the previous values. FLINT runs it on integers over the common denominator (n + 1)\#, so that no GCD is needed until the final division by the true denominators.

• Power series modulo primes (_arith_bernoulli_number_vec_multi_mod, n \ge 700). The even Bernoulli numbers are the coefficients of

  \frac{x^2}{\cosh x - 1} = \sum_{k \ge 0} \frac{(2 - 4k)\, B_{2k}}{(2k)!} x^{2k},

  obtained by one power series inversion of (\cosh x - 1)/x^2, which has only the even terms [BuhlerCrandallSompolski1992]. Multiplying by the known denominators and factorials modulo each prime gives the numerators modulo the prime, and Chinese remaindering the numerators.

Euler numbers#

The Euler numbers E_n, 1/\cosh t = \sum E_n t^n/n!, are integers with |E_{2m}| = 2^{2m+2} (2m)!\, \beta(2m+1)/\pi^{2m+1}, where \beta(s) = \sum_{k \ge 0} (-1)^k (2k+1)^{-s} is the Dirichlet beta function. FLINT computes an isolated E_n from this formula for n < 2000 and with a multimodular method combined with the beta function above (arb_fmpz_euler_number_ui), and a vector by inverting the power series of \cosh modulo primes (arith_euler_number_vec).

Use in FLINT#

The generic functions call these: gr_bernoulli_ui uses bernoulli_fmpq_ui over \mathbb{Q} (and in rings without a real precision, followed by a conversion), gr_bernoulli_vec uses bernoulli_fmpq_vec_no_cache over \mathbb{Q}, and gr_eulernum_vec uses arith_euler_number_vec over \mathbb{Z}; in real and complex rings they compute numerically with arb.