Verifying an Efficient Algorithm for Computing Bernoulli Numbers

For presentation purposes, we start by describing the algorithm for the remainder, which works for arbitrary $n$. The algorithm for chunk-aligned $n$ uses the same implementation techniques with the addition of a tabulation optimisation, which we describe in Section 3.2.

For the aligned part of the summation, we apply another crucial optimisation that is also described by Harvey. Instead of iterating through the sum in single steps, we iterate in chunks of size $w{w}^{\prime }$. Note that $\sigma (p,2^{j}s)$ being $1$ or $-1$ is equivalent to the $j$th digit of the binary expansion of $s/p$ being $0$ or $1$, which allows us to determine the values of $\sigma (p,2^{j}s)$ for $w{w}^{\prime }$ consecutive values of $j$ in a single step by computing the corresponding $w{w}^{\prime }$ bits of the binary expansion of $s/p$. For this we define

and rewrite the inner sum for multiples of $w{w}^{\prime }$ as

where $\cdot !k$ selects the $k$th bit of its argument, such that $\text{blbits}(i,j)!k=\text{digf}(2,iw{w}^{\prime }+jw+k)$. The table is defined as

Since $n$ is typically much larger than $w$ and $2^w$, most of the computational work goes into computing the table, which our optimisations make cheap: in an outer loop, we obtain chunks of $w{w}^{\prime }$ bits of the binary expansion and incrementally maintain $x^{iw{w}^{\prime }}$. The inner loop only bit-shifts and masks these bits, indexes into the table, and adds $x^{iw{w}^{\prime }}$ to table entries. Additionally, the implementation makes similar optimisations to those we described in Sec. 3. In particular, we use the relaxed Montgomery form for the table entries to allow for a mostly branch free implementation of the inner loop.⁵⁵5The computation of the binary expansion of $s/p$ we describe below does use a branch to correct a possible ‘off-by-one’ error. However, this happens so rarely that due to branch prediction, the cost of this branch is negligible. Moreover, the size of the table ($w^{\prime }{2}^w$ double-width machine words, i.e. $16\mathrm{KiB}$ in our current setting) is chosen such that it fits comfortably inside the L1d cache of a typical modern CPU.

To incrementally obtain the bits of the binary expansion, we maintain the value of modf, and observe that $\text{modf}(b,x,p,i+1)=(b\cdot \text{modf}(b,x,p,i))modp$. However, instead of performing an expensive division by $p$, we apply a further optimisation: we define $p^{-}=\lfloor 2^{128}/p\rfloor$. Then, for and , we have $\lfloor 2^{64}x/p\rfloor \in \{\lfloor p^{-}x/2^{64}\rfloor ,\lfloor p^{-}x/2^{64}\rfloor +1\}$. Using this, we can reduce the division and modulo operation to a multiplication operation. To be able to detect the $+1$ case, we limit , and define:

Here, >> denotes the left shift operation and AND the bitwise ‘and’. With this, we get:

This shows how to obtain the next digit and maintain the value of modf. Note that $\lfloor 2^{128}/p\rfloor$ is precomputed outside the loop. Figure 2 displays the algorithm for computing the table.

After some initialisation, the outer loop iterates over the chunks. In each iteration of the outer loop, the binary expansion of $s/p$ is maintained and another 64 bits of digits are obtained in cbits. The inner loop iterates over the blocks, obtains the current block bits, and updates the table. After the inner loop, the outer loop maintains the value curx of $x^{iw{w}^{\prime }}$.

Finally, we assemble the aligned inner sum, the computation of the remainder, and the fallback to the slower algorithm when $2^k{\equiv }_{p}0$ (i.e. when $kmod{\text{ord}}_p(2)=0$) to obtain the algorithm harvey_select2 with the following correctness theorem:

That is, under the listed preconditions, harvey_select2 will return the Montgomery form of an integer $r$ with $r{\equiv }_p{B}_k$.

Next, we address some of the preconditions. First, we compute $p^{\prime }=p^{-1}mod{2}^W$ via Hensel lifting. Next, as mentioned before, when $(p-1)\mid k$ we have $N_k{\equiv }_p-D_k/p$, which allows us to easily determine $N_{k}modp$ directly. Moreover, when $k\ge p-1$, we apply the previously-mentioned range reduction based on Kummer’s congruence. Finally, we multiply the result with $D_k$ to obtain $N_{k}modp$, and convert the result from Montgomery form to normal numbers.

The resulting algorithm is displayed in Fig. 3.

In addition to the parameters $p$, $g$, $\mathit{op2}$, and $k$, it also takes the prime factorisation $\mathrm{𝑑𝑓𝑠}$ of $D_k$. The function prod_list_exc2 computes $D_k/p$ by multiplying all prime factors except $p$. The correctness theorem is:

In words: the algorithm returns $N_{k}modp$ for any even natural number $k\ne 0$ and any odd prime less than $2^{31}$.

The function estimate_check_a_priori_bound k (cf. Fig. 4) returns a pair $(ok,Y+1)$. If $ok=\text{True}$, then $Y$ is an upper bound for the largest prime number that will be needed to reconstruct $N_k$. If $ok=\text{False}$, we cannot guarantee that our 32-bit arithmetic will not overflow, and we abort the algorithm. The first ingredient to compute $Y$ is a bound on $B_k$:

The prime sieving algorithm mk_factor_prime_cache K returns a tuple $(\mathrm{𝑓𝑐},\mathrm{𝑝𝑐})$, where $\mathrm{𝑝𝑐}$ (“prime cache”) is the list of primes $p$ with , and $\mathrm{𝑓𝑐}$ (“factor cache”) encodes a mapping that maps any integer between $0$ and $K$ to its smallest prime factor (if it is composite) or to 0 (if it is prime or $\le \mathrm{\hspace{0.17em}1}$). The sieving algorithm starts by mapping all numbers to $0$ and then proceeds in the usual fashion. When sieving is finished, the list of primes is extracted from the factor cache.

The advantage of the factor cache is that it allows us to not only quickly list prime numbers and determine whether a number is prime, but also to factor numbers efficiently. This speeds up the compute_harvey_pgos algorithm (cf. Sec 4.4).

We can now determine the prime factorisation of $D_k$ by simply taking all primes $p\le k$ with $(p-1)\mid k$. The algorithm adjust_primes2 dfs k pc sums up ${\log }_{2}p$ for all of these to get an approximation of ${\log }_2{D}_k$ that easily fits in a machine integer. Together with our estimate for ${\log }_2|B_k|$, we can now compute a relatively precise estimate of ${\log }_2|N_k|$.

We then choose a prefix $P$ of the sieved primes such that ${\sum }_{p\in P}{\log }_2(p)\ge {\log }_2|N_k|$. For this, we use a fixed-point approximation of ${\log }_{2}p$. This avoids using floating-point numbers or arbitrary precision integers and is still precise enough in practice. Note that we do not actually store the prime 2 in the sieve or pruned list of primes, but still consider it for the bound. Before reconstructing $|N_k|$ (cf. Sec. 4.6), we do add the congruence $N_k{\equiv }_{2}1$ (which holds for any even $k$) to our modular information.

The algorithm compute_harvey_pgos pc fc takes the pruned prime list and the factor map and computes a list of task entries of the form $(p,g,o,0)$. Here, $p$ is the prime itself, $g$ is a generator of $\mathbb{Z}/p\mathbb{Z}$, $o$ is the multiplicative order of 2 in $\mathbb{Z}/p\mathbb{Z}$, and $0$ is a placeholder that will be filled in by Harvey’s algorithm.

A number $g$ is a generator of $\mathbb{Z}/p\mathbb{Z}$ iff it has maximal order, i.e. if ${\text{ord}}_p(g)=p-1$. We can determine whether $g$ is a generator by checking that $g^{(p-1)/q}{\not\equiv }_{p}1$ for all prime factors $q$ of $p-1$. We find the smallest generator $g_p$ by simply checking $g=2,3,\mathrm{\dots },p-1$ and taking the first one that works.⁶⁶6Empirically, $g_p$ is almost always quite small. For the first ${10}^6$ primes, the mean and maximum of $g_p$ are 4.91 and 94, with $g_p\in \{2,3\}$ in 60 % of the cases. It is known that $g_p\in O(p^{\frac{1}{4}+\epsilon })$ and the generalised Riemann hypothesis implies $g_p\in O({\log }^{6}p)$.

As for computing the order: ${\text{ord}}_p(2)$ is defined as the smallest positive integer $n$ such that $2^n{\equiv }_{p}1$. Lagrange’s theorem tells us that ${\text{ord}}_p(2)\mid p-1$. Therefore, every factor $q$ in the prime factorisation of ${\text{ord}}_p(2)$ must also be a factor of $p-1$, and we can determine its multiplicity ${\nu }_q({\mathrm{ord}}_p(2))$ by letting $i\leftarrow 0$, $x\leftarrow 2^{(p-1)/q^e}modp$ (where $e={\nu }_q(p-1)$) and then repeating $i\leftarrow i+1$, $x\leftarrow x^{q}modp$ until $x=1$, at which point $i={\nu }_q({\mathrm{ord}}_p(2))$ as desired.

We are not aware of better algorithms for computing the order and finding a generator; however, even the worst-case running time of these is negligible compared to the main part of our algorithm.

To compute the modular information in parallel, we instantiate a generic parallel map combinator with bernoulli_num_harvey_wrapper2 (cf. Sec. 3.3). We elide the boilerplate code for the instantiation, and only display the resulting correctness lemma in Fig. 5.

Figure 6 displays our Isabelle formalisation of Harvey’s full algorithm. After handling some special cases for $k=0$, $k=1$, and odd $k$, it computes the modular information and a prime factorisation of the denominator. If the $\mathrm{𝑜𝑘}$ flag is false, we abort the algorithm. Otherwise, the $(p,g,o,r)$ tuples in $\mathrm{𝑝𝑔𝑜𝑟𝑠}$ are converted to $(r,p)$ pairs in $\mathrm{𝑎𝑚𝑠}$. We also add the pair $(1,2)$ since $N_k$ is always odd. At this point, we know that for each pair $(r,p)\in \mathrm{𝑎𝑚𝑠}$, we have $N_{k}modp=r$. Moreover, the primes $P$ in these pairs are disjoint, and their product is greater than $N_k$. Chinese remaindering gives us a pair $(\mathrm{𝑛𝑢𝑚},m)$ with $\mathrm{𝑛𝑢𝑚}=N_{k}modm$ and $m={\prod }_{p\in P}p$. Thus, $N_k=\mathrm{𝑛𝑢𝑚}$ if $N_k\ge 0$ and $N_k=\mathrm{𝑛𝑢𝑚}-m$ if , and we know that iff $4\mid k$. Finally, we compute the denominator $\mathrm{𝑑𝑒𝑛𝑜𝑚}=D_k$ from its prime factorisation.