A multimodular algorithm for computing Bernoulli numbers

The Bernoulli numbers $B_{0},B_{1},B_{2},\ldots$ are rational numbers defined by

In this paper we focus on the problem of computing $B_{k}$, as an exact rational number, for a single large (even) value of $k$. To date, the most efficient algorithms for computing an isolated $B_{k}$ have been based on the formula

where $\zeta(s)$ is the Riemann zeta function [IR90, p. 231]. One first computes a high-precision approximation for $B_{k}$ via (1), using the Euler product for $\zeta(s)$. The exact denominator of $B_{k}$ is known by the von Staudt–Clausen theorem. Provided that the numerical approximation to $B_{k}$ has enough correct bits — it turns out that $O(k\log k)$ bits is enough — we may recover the numerator of $B_{k}$ exactly. We will refer to this algorithm and its variants as the ‘zeta-function algorithm’.

The zeta-function algorithm has been rediscovered numerous times. The formula (1) itself goes back to Euler. Chowla and Hartung [CH72] give an ‘exact formula’ for $B_{k}$, the evaluation of which is tantamount to executing the algorithm, but using the defining sum for $\zeta(k)$ instead of the Euler product, and using $2(2^{k}-1)$ in place of the exact denominator of $B_{k}$. Fillebrown [Fil92] mentions [CH72] and points out that using the Euler product is asymptotically faster; she gives a detailed time and space complexity analysis, and attributes the algorithm to Herbert Wilf (unpublished). Fee and Plouffe [FP07] state that they discovered and implemented the zeta-function algorithm in 1996. Kellner [Kel02] also refers to [CH72] and suggests using the Euler product. He mentions some data produced by Fee and Plouffe, although not their algorithm. According to Pavlyk [Pav08], the implementation in Mathematica 6 [Wol07] derives from Kellner’s work. Stein [Ste07, p. 30] indicates another independent discovery by Henri Cohen, Karim Belabas, and Bill Daly (unpublished), that led to the implementation in PARI/GP [Par07].

In this paper we present a new algorithm for computing $B_{k}$ that does not depend on (1) or any properties of $\zeta(s)$. The basic idea is to compute $B_{k}$ modulo $p$ for sufficiently many small primes $p$, and then reconstruct $B_{k}$ using the Chinese Remainder Theorem. Using a parallel implementation, we have computed $B_{k}$ for $k=10^{8}$, improving substantially on the previous record of $k=10^{7}$ [Pav08].

Let $M(n)$ denote the time (bit operations) required to multiply two $n$-bit integers. We may assume that $M(n)=O(n\log^{1+\varepsilon}n)$, using for example the Schönhage–Strassen algorithm [SS71]. According to Fillebrown’s analysis [Fil92, p. 441], the zeta-function algorithm for computing $B_{k}$ requires $O(k)$ multiplications of integers with $O(k\log k)$ bits, so its running time is $O(k^{2}\log^{2+\varepsilon}k)$. The new algorithm also runs in time $O(k^{2}\log^{2+\varepsilon}k)$. However, we will see that over the feasible range of computation the running time behaves more like $O(k^{2}\log k)$.

The new algorithm is easily parallelisable, as the computations for the various $p$ are independent. Only the final stages of the modular reconstruction are more difficult to execute in parallel, but these account for only $O(k^{1+\varepsilon})$ of the running time. During the modular computations, each thread requires only $O(\log^{1+\varepsilon}k)$ space. The zeta-function algorithm may also be parallelised, for instance by mapping Euler factors to threads, but this requires $O(k\log k)$ space per thread.

Throughout this section $p\geq 5$ denotes a prime, and $k$ an even integer satisfying $2\leq k\leq p-3$. Our algorithm for computing $B_{k}$ modulo $p$ is based on the following well-known congruence. Let $0<c<p$, and suppose that $c^{k}\not\equiv 1\pmod{p}$. Then

where

Equation (2) may be deduced by reading the statement of Theorem 2.3 of [Lan90, §2] modulo $p$ (as is done in the proof of Corollary 2 to that theorem).

Equation (2) may be used to compute $B_{k}$ modulo $p$ directly, but computing $x^{k-1}$ modulo $p$ for each $x$ is unnecessarily expensive. Instead, we select a generator $g$ of the multiplicative group $(\mathbf{Z}/p\mathbf{Z})^{\times}$, put $c=g$, and rewrite the sum as

Note that $g^{k}\not\equiv 1\pmod{p}$ since $p-1\mathbin{\nmid}k$. For $x\not\equiv 0\pmod{p}$ we have $h_{c}(-x)=-h_{c}(x)$, and $g^{(p-1)/2}\equiv-1\pmod{p}$, so we obtain the half-length sum

leading to the following algorithm.

In this section we present a multimodular algorithm for computing $B_{k}=N_{k}/D_{k}$ as an exact rational number. Before getting to the algorithm proper, we make some preliminary calculations. In what follows, we assume that $k\geq 4$ and that $k$ is even.

The denominator $D_{k}$ is given by the von Staudt–Clausen theorem [IR90, p. 233]:

Note that $D_{k}\leq 2^{k+1}$, since $D_{k}$ divides $2(2^{k}-1)$ [CH72, p. 114].

The size of the numerator $N_{k}$ may be bounded quite tightly as follows. From (1) and Stirling’s approximation [Lan97, p. 120] we have

Since $\zeta(k)\leq\zeta(4)$ and $12k\geq 48$, a short calculation shows that $|N_{k}|<2^{\beta}$ where $\beta:=\lceil(k+0.5)\log_{2}k-4.094k+2.470+\log_{2}D_{k}\rceil$.

For any $X>0$, let $M_{X}:=\prod_{p\leq X,p-1\mathbin{\nmid}k}p$. Our strategy will be to select $X$ so that $M_{X}\geq 2^{\beta}$, and then compute $N_{k}$ modulo $M_{X}$. This ensures that we have sufficient information to reconstruct $N_{k}$. In the algorithm itself our choice of $X$ will be quite sharp, but we also need a rough a priori estimate. We will take $Y:=\max(37,\lceil(k+0.5)\log_{2}k\rceil)$. To check that this works, we will use the estimate $\sum_{p\leq x}\log p\geq 0.73x$, a weak form of the prime number theorem, valid for $x\geq 37$ [Nar00, p. 111]. Then we have

so that $M_{Y}\geq 2^{\beta+2}$. Algorithm 2 presents the resulting algorithm; several important details are given in the proof of the following theorem.

In this section we sketch several techniques that yield a large constant factor improvement in the speed of computing $B_{k}$ modulo $p$, for almost all $k$ and $p$. In experiments we have found that these techniques yield an increase in speed by a factor in excess of $150$ compared to an efficient implementation of Algorithm 1.

We first perform some preparatory algebra. Consider again the congruence (2). Let $n$ be the multiplicative order of $c$ in $(\mathbf{Z}/p\mathbf{Z})^{\times}$, and put $m=(p-1)/n$. Let $g$ be a generator of $(\mathbf{Z}/p\mathbf{Z})^{\times}$. Then $\{g^{i}c^{-j}\}_{0\leq i<m,0\leq j<n}$ forms a complete set of representatives for $(\mathbf{Z}/p\mathbf{Z})^{\times}$. Putting $x=g^{i}c^{-j}$, we obtain

Considerations similar to those of §2 allow us to halve the number of terms, as follows. First suppose that $n$ is even. Since $h_{c}(-x)=-h_{c}(x)$ and $c^{n/2}\equiv-1\pmod{p}$, the double sum in (5) becomes

Now suppose that $n$ is odd. Then $m$ is even, and we have $(-g^{m/2})^{n}\equiv(-1)^{n}g^{(p-1)/2}\equiv(-1)(-1)\equiv 1\pmod{p}$, so that $-g^{m/2}\equiv c^{\lambda}\pmod{p}$ for some $0\leq\lambda<n$. The double sum in (5) becomes

We combine both cases into a single statement by writing

where

The sum in (6) may be evaluated by an algorithm similar to Algorithm 1, with an outer loop over $i$ and an inner loop over $j$; we omit the details. The inner loop is executed $n^{\prime}m^{\prime}=(p-1)/2$ times, the same number of iterations as in the main loop of Algorithm 1. Note that if $c$ is chosen ‘randomly’, then we expect $n^{\prime}$ to be quite large, so it is imperative to make the inner loop as efficient as possible.

To this end, we specialise to the case $c=1/2\pmod{p}$; this will allow us to exploit the fact that multiplication by $2$ modulo $p$ is very cheap. Of course, this specialisation does not work if $2^{k}\equiv 1\pmod{p}$. In practice, such primes tend to be rare for any given $k$, and we may fall back on Algorithm 1 to handle them.

Directly from the definition of $h_{c}(x)$, we have $h_{1/2}(x)=-f(x)/4$ where

From (6) we obtain

where $n$ is the order of $2$ in $(\mathbf{Z}/p\mathbf{Z})^{\times}$, and $n^{\prime}$ and $m^{\prime}$ are defined as before.

In the algorithm that evaluates (7), the inner loop is considerably simpler than the main loop of Algorithm 1. Given $g^{i}2^{j}$ at the beginning of the loop, we compute the next term $g^{i}2^{j+1}$ by doubling it and subtracting $p$ if necessary, and $f(g^{i}2^{j})$ is $-1$ or $+1$ according to whether the subtraction occurred. Therefore only a single modular multiplication is required on each iteration, to keep track of $(2^{k-1})^{j}$.

Next we describe further optimisations for computing a sum of the form

(For evaluating (7), we take $N=n^{\prime}$, $r=2^{k-1}$, $s=g^{i}$.) The following observation is crucial: if $0<s<p$, and if the binary expansion of $s/p$ is $0.b_{0}b_{1}b_{2}\ldots$ (where each $b_{j}$ is 0 or 1), then $f(2^{j}s)=(-1)^{b_{j}}$. Indeed, the binary expansion of $2^{j}s/p$ is obtained by shifting that of $s/p$ to the left by $j$ digits, so

and then

Assume that we have available a routine for computing the binary expansion of $s/p$, whose output is a sequence of base-$2^{m}$ digits. Strip off the last $N\bmod m$ terms of (8) and compute them separately; we may then assume that $N$ is divisible by $m$. Put $j=mt+u$, so that (8) becomes

We may now use a caching strategy as follows. Maintain a table $\{T_{\sigma}\}$ of length $2^{m}$, where $\sigma\in\{+1,-1\}^{m}$. Each $T_{\sigma}$ is a residue modulo $p$, and is initialised to zero. For each $0\leq t<N/m$, add $(r^{m})^{t}$ to the table element $T_{\sigma}$, where $\sigma$ is determined from the $t$-th base-$2^{m}$ digit of $s/p$ by

(In practice, $\sigma$ is represented by the sequence of bits themselves, and is obtained by a simple bit-mask.) After finishing the loop over $t$, compute the $2^{m}$ values $V_{\sigma}=\sum_{0\leq u<m}r^{u}\sigma_{u}$, where $\sigma=(\sigma_{0},\ldots,\sigma_{m-1})$. Then the desired sum is given by $\sum_{\sigma}V_{\sigma}T_{\sigma}$. The main benefit of this approach is that in the inner loop we only perform $N/m$ modular multiplications, instead of $N$. The tradeoff is that we must maintain the table $\{T_{\sigma}\}$, and some extra work is required at the end to assemble the information from the table. For this to be worthwhile, we should have $N\gg 2^{m}$. We should also choose $m$ so that the base-$2^{m}$ digits of $s/p$ can be extracted efficiently. Moreover, memory locality problems can offset the savings in arithmetic if $m$ is too large. We have found in practice that $m=8$ works quite well.

In the implementation discussed in §5, we made several further optimisations, which we describe only briefly.

For better memory locality — and indeed less total memory consumption — we split (9) into blocks and process each block separately. This avoids needing to store and then later retrieve too many bits of the expansion of $s/p$. Instead of computing $s/p$ for each $s$ separately, we precompute an approximation to $1/p$, and then multiply it by each $s$ that occurs, taking advantage of the speed of multiplication by a single-word integer compared to the corresponding division.

We use several tables instead of just one. For example, on a 32-bit machine, we use four tables: the highest 8 bits of each word of $s/p$ contribute to the first table, the next 8 bits to the second table, and so on. This reduces further the number of modular multiplications. When adding values into a table, we skip the reduction modulo $p$, delaying this until the end of the loop. Of course, this is only possible when the word size is large enough compared to $p$, otherwise overflow may occur. Finally, instead of re-initialising the table on each iteration of the outer loop, we simply continue to accumulate data into the same table. This partly mitigates against the overhead incurred when $m^{\prime}$ is unusually large compared to $n^{\prime}$.

Most of the modular arithmetic uses Montgomery’s reduction algorithm [Mon85].

The author implemented the above algorithms in a C++ package named bernmm (Bernoulli multi-modular). The source code is released under the GNU General Public License (GPL), and is freely available from the author’s web page. It depends on the GMP library [Gra08] for arbitrary precision arithmetic, and on NTL [Sho07] for some of the single-precision modular arithmetic. It is included as a standard component of the Sage computer algebra system (version 3.1.2) [SJ05], accessible via the bernoulli function.

The parallel version distributes the modular computations among the requested number of threads, and depends on the standard pthreads API. For the reconstruction phase, the lowest layers of the subproduct tree are performed in parallel, but the number of threads decreases towards the top of the tree, since the extended GCD routine is not threaded.

Table 1 shows timing data for $10^{4}\leq k\leq 10^{8}$, spaced by intervals of $\sqrt{10}$, for bernmm and several existing implementations of the zeta-function algorithm. The timings were obtained on a 16-core 2.6GHz AMD Opteron (64-bit) machine with 96 GB RAM, running Ubuntu Linux. The last two rows of the table ($k=\lfloor 10^{7.5}\rfloor$ and $k=10^{8}$) improve on the prior record $k=10^{7}$ set by Pavlyk [Pav08]. The values $B_{31622776}$ and $B_{10^{8}}$ may be downloaded from the author’s web page.

The compiler used was gcc 4.1.3. We used the version of GMP shipped with Sage 3.1.2, which is the same as GMP 4.2.2, but also includes assembly code written by Pierrick Gaudry that improves performance on the Opteron, and code by Niels Möller that implements a quasi-linear time extended GCD algorithm [Möl08]. We used NTL version 5.4.1, also included in Sage.

The timings for PARI/GP were obtained for version 2.3.3, using the same patched version of GMP. We used the bernfrac function from the GP interpreter. The timings for calcbn, a specialised program for computing Bernoulli numbers, were obtained from calcbn 2.0 [Kel08], again using the same patched GMP. The timings for Mathematica were obtained from the Linux x86 (64-bit) version of Mathematica 6.0, using the BernoulliB function. Neither PARI/GP nor Mathematica supply a parallel implementation as far as we are aware.

The peak memory usage for computing $B_{10^{8}}$ was approximately 5 GB; this occurred during the final extended GCD operation. By comparison, PARI/GP used 4.3 GB for $k=10^{7}$, and calcbn used 1.7 GB for $k=10^{7}$ (with ten threads).

For each implementation, we observe that the ratio of the times in adjacent rows of the table is approximately 10, reflecting the predicted quasi-quadratic asymptotic running time of both the zeta-function algorithm and the multimodular algorithm.

We checked our results by two methods. First, we checked that as a real number, the proposed value for $B_{k}$ agrees with the estimate $|B_{k}|\approx 2(2\pi)^{-k}k!$. Since we computed $B_{k}$ by modular information alone, this is a very strong consistency check. Second, we reduced modulo $p$ the proposed value for $B_{k}$, and compared this against the output of Algorithm 1, for a few primes distinct from those primes that we used to compute $B_{k}$. Again, this is a very strong consistency check — stronger in fact, since it involves all of the bits of the proposed $B_{k}$.

Mathematica, PARI/GP and calcbn use GMP internally, so naturally any improvement in GMP’s large integer arithmetic will improve the timings reported in the table. In this connection we mention the improvements in multiplication speed reported by [GKZ07]; their paper also gives comparisons with other implementations of large-integer arithmetic.

Many thanks to Joe Buhler, Bernd Kellner, and Andrew Sutherland for their comments on a draft of this paper, and to the Department of Mathematics at Harvard University for providing the hardware on which the computations were performed.