A multimodular algorithm for computing Bernoulli numbers

Math. Comp. 79 (2010), 2361–2370 (DOI).

arXiv preprint (October 2008).

Abstract

We describe an algorithm for computing Bernoulli numbers. Using a parallel implementation, we have computed B_k for k = 10⁸, a new record. Our method is to compute B_k modulo p for many small primes p, and then reconstruct B_k via the Chinese Remainder Theorem. The asymptotic time complexity is O(k² log^2+ε k), matching that of existing algorithms that exploit the relationship between B_k and the Riemann zeta function. Our implementation is significantly faster than several existing implementations of the zeta-function method.

Implementation

I implemented the algorithm described in the paper in the bernmm package.

Data

The 10,000,000-th Bernoulli number

• The denominator of B_10,000,000 is 9601480183016524970884020224910.
• Download the rational number B_10,000,000 as a Sage object (24MB).
• Download the numerator of B_10,000,000 in decimal as plain text (24MB).
• The value may be verified against Pavlyk's computation (see also this blog post).

The 31,622,776-th Bernoulli number

• The denominator of B_31,622,776 is 369037807590.
• Download the rational number B_31,622,776 as a Sage object (83MB).
• Download the numerator of B_31,622,776 in decimal as plain text (81MB).

The 100,000,000-th Bernoulli number

• The denominator of B_100,000,000 is 394815332706046542049668428841497001870.
• Download the rational number B_100,000,000 as a Sage object (284MB).
• Download the numerator of B_100,000,000 in decimal as plain text (278MB).