"Practical" Sieve of Eratosthenes from "Primes Numbers - A Computational Perspective"

Question

Consider the following pseudocode for the Sieve of Eratosthenes, giving us the primes up to $N$:

1) List the numbers $2$ to $N$.

2) Let $p=2$.

2) Cross out $p^2$, then cross out $(p+1)p, (p+2)p, (p+3)p,..., (p+j)p$, where $j$ is the smallest integer such that $(p+j+1)p > N$.

3) Set $p$ equal to the first number in the list which isn't crossed out. If $p^2 > N$, then we stop, otherwise, go to step 2)

If you follow the above algorithm, then the numbers which are crossed out are the composite numbers between $2$ and $N$, and the remaining numbers are prime. I understand why this works, and why it's reasonably fast (because step 2 can be achieved only with one squaring and $j$ additions). This is pretty much the standard sieve of Eratosthenes, which you can find in various textbooks, as well as Wikipedia and Wolfram Mathworld.

However, I've recently started reading "Prime Numbers - A Computational Perspective" by Crandall and Pomerance. This is an excellent textbook, but I'm confused about Algorithm 3.2.1, which they call the "Practical Eratosthenes sieve" (page 122 in the second edition). I can see the resemblance to the algorithm above, but I don't understand exactly what's going on, why certain restrictions are assumed, and why it works. I've copied the algorithm (almost) verbatim here (I've clarified some of the ambiguous formatting). It is written in a fairly standard, C-style syntax:

"Algorithm 3.2.1 (Practical Eratosthenes sieve).

This algorithm finds all primes in an interval $(L, R)$ by establishing Boolean primality bits for successive runs of $B$ odd numbers. We assume $L$, $R$ even, with $R>L$, $B | R − L$ and $L> P = \lceil{\sqrt{R}}\rceil$. We also assume the availability of a table of the $\pi(P)$ primes $pk ≤ P$.

1) [Initialize the offsets]

for($k \in [2, \pi(P)]$) {

$\quad q_k = (-\frac{1}{2}(L +1+ p_k)) \bmod p_k$;

}

2) [Process blocks]

$T = L$;

while($T <R$) {

$\quad$for($j \in [0, B − 1]$) {

$\quad\quad b_j = 1$;

$\quad$}

$\quad$for($k \in [2, \pi(P)]$) {

$\quad\quad$for($j = q_k$; $j<B$; $j = j + p_k$) {

$\quad\quad\quad b_j = 0$;

$\quad\quad$}

$\quad\quad q_k = (q_k − B) \bmod p_k$;

$\quad$}

$\quad$for($j \in [0, B − 1]$) {

$\quad\quad$if($b_j$ == 1) report $T + 2j + 1$; // Output the prime $p = T + 2j + 1$.

$\quad$}

$\quad T = T + 2B$;

}"

(I apologise for my terrible formatting. I couldn't think of a better way of doing it.)

1) What is the significance of $B$? Is it just to lower memory requirements? Is there a "sweet spot" for choosing a value of $B$? Obviously, we could choose $B=1$ or $B = R-L$, but what should we consider when choosing a value of $B$? Additionally, why do we need $B | R - L$?

2) Why do we need $L > \lceil{\sqrt{R}}\rceil$?

3) What are the $q_k$ for? Why is their initial value set to $q_k = -(L + 1 + p_k)/2 \bmod p_k$?

4) What is happening in the main while loop? I can see that by setting the $b_j = 0$, we're "crossing out" composites as before, but I don't understand any of the logic (why increment $j$ by $p_k$ each time? Why set $q_k = (q_k - B) \bmod p_k$? Why set $T = T + 2B$?).

I'm sorry if this is a terribly formed question, but I don't even know where to begin with this. To me, it only slightly resembles the naive sieve of Eratosthenes mentioned at the start of my post, and I don't see why the Algorithm 3.2.1 should be more efficient or more appropriate.

I've asked a lot, but if anyone can tackle any of my queries, it'd be massively appreciated.

Answer 1

1) I can't entirely speak for the authors, but yes the block size $B$ looks to be chosen primarily by memory considerations. Modern computer systems have three or more levels of working memory which trade off capacity for latency: persistent storage (hard drive/SSD), RAM, and on-board cache (some of these levels are themselves split into more sublevels: multiple levels of cache or hybrid SSD/HD drives, each affording more trade-offs).

Abstractly, choosing your working set (which will be proportionate to $B$ because the algorithm makes a lot of writes all across vector $b$ rather than spending time in any local section) so that it fits entirely into one level of memory will grant it the speed benefits of that tier, which is roughly an order of magnitude faster than the next more-capacious tier. An algorithm that fits into L1 cache can perform basic read/write operations vastly faster than one that frequently swaps to hard disk. And sure, if $B$ is too large to fit in external storage, then it won't run at all.

On the other hand, making $B$ too small incurs the overhead of scanning through the list of primes $(R-L)/B$ times, so there is strong benefit in making $B$ as large as possible while fitting into a given tier (this is what makes the Sieve much faster than testing many individual numbers).

Since there are two opposing forces pushing $B$ towards small and large values, a sweet spot lies somewhere between the extremes. Note also that the cost of making $B$ small depends on the size of $R$, so the sweet spot depends on the size of the problem and not entirely on the hardware (which determines the cost of making $B$ large).

Making $B \mid R-L$ seems like just a minor convenience, to avoid handling partial blocks. The algorithm as written would be incorrect for $B \nmid R-L$ since the final iteration checks values $>R$ against a table of primes that only goes up to $\sqrt{R}$, although that is easy to account for. In fact I believe this is a typo, and it should be $B\mid (R-L)/2$.

2) We choose $P = \lceil \sqrt{R} \rceil$ for the reason just described, and we choose $L>P$ to simplify the logic: a number $n \in (L,R)$ is prime iff no prime $p_k$ divides $n$. If $L\le P$ we'd have to amend this with "except for when $n = p_k$".

I believe the rest of your questions are best answered by describing the algorithm at a high-level, which should indirectly answer most of them simultaneously:

1. The idea is to go partition the interval $(L,R)$ into blocks of size $2B$, of the form $(L+2Bm,L+2B(m+1))$, and sieve each block one at a time for $m=0, 1, 2, \ldots$. The variable $T$ encodes the current value of $L+2Bm$ (hence $T \leftarrow T+2B$).

2. For each block, we only need to consider the odd values, so there are only $B$ numbers to keep track of instead of $2B-1$. Since we are scanning over odd numbers, we are really crossing off values that are $2p_k$ apart rather than $p_k$, but since the bit vector $b$ corresponds only to odd numbers, the spacing within $b$ is still $p_k$.

3. Since our block starts at $T$ rather than $0$, we need to calculate where in the block we can find the first odd multiple of $p_k$ (so we know where to start crossing off numbers), and that's the purpose of $q_k$. For the first block, $T=L$ so the first block represents the $B$ odd values $L+1, L+3, L+5, \ldots, L+2B-1$. The initial value of $q_k$ represents the index of the smallest multiple of $p_k$ in that list (try working it out yourself so you get a feel for the index arithmetic and where the factor of $\tfrac12$ arises, why it depends on $p_k$, and where you end up using the assumption that $L$ is even).

4. When we move to the second block, we are now considering $L+2B+1, L+2B+3, \ldots, L+4B-1$, and the index of the first multiple of $p_k$ shifts by a predictable amount. This explains the $q_k \leftarrow (q_k-B)\text{ mod } p_k$. The same shift occurs happens for every subsequent block.