A multimodular algorithm for computing Bernoulli numbers

Abstract

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