The sieve of Eratosthenes

To issue a 2048-bit RSA key, a program needs two random 1024-bit primes, each with 309 decimal digits. There is no ready list of such numbers and there cannot be, so the program takes a random odd number and goes up from it until it meets a prime. Among odd numbers of this size a prime occurs on average one in 355, and checking a candidate with the Miller–Rabin test means raising to a power with a 309-digit exponent modulo a 309-digit number, more than a thousand multiplications of huge numbers. You can check them all in a row. Or you can sieve first: cross out the candidates that have a small divisor (more than eighty per cent of them) and spend the expensive test only on the survivors. Below this happens in your browser, on real numbers.

Each cell is an odd candidate: x, x + 2, x + 4 and so on. First the sieve crosses out every third cell, the multiples of three, then every fifth, seventh, and after a few primes the window is already full of holes. Then the expensive test goes over the survivors, and the first number it does not refute is prime. The test is probabilistic, but after 64 rounds the chance of error is below 4⁻⁶⁴. Move the “sieve up to” slider. Without a sieve every candidate gets the test. A sieve up to a thousand leaves about 16% of candidates, and about six times fewer expensive checks are needed. A sieve up to a million leaves 8%: the gain grows, but ever more slowly, while the sieve itself gets noticeably more expensive: it is 78 thousand divisions of a huge number by small ones.

The idea

Do not test numbers for divisibility one at a time; go the other way, from divisors to multiples. Each prime p shows by itself where its multiples stand: every p numbers. Stepping through the array with step p and making marks takes not a single division, only additions.

The classic sieve finds all primes up to n. Write out the numbers from 2 to n. The smallest number not crossed out is prime: if it had a smaller divisor, it would have been crossed out as a multiple of that divisor. Declare it prime, cross out its multiples and repeat. Two details make the sieve fast:

• the multiples of p are crossed out starting from p², not from 2p: the numbers 2p, 3p, …, (p − 1)p have a divisor smaller than p and are already crossed out;
• stop when p² exceeds n: there is nothing left to cross out, and everything remaining is prime.

A detailed school-level walkthrough with a picture and a proof is in the chapter “Prime numbers” of “Mathematics, the Queen of the Sciences”. Here we are interested in something else: how much work the sieve does, where it does it for nothing, and how it is sped up in practice.

How many times one number is crossed out

Two crosses out 4, 6, 8, …, three crosses out 9, 12, 15, … The number 12 falls under both: it is crossed out twice. The number 30, three times: by two, three and five. The sieve does not remember that a number is already crossed out and wastes time on repeats. Go through the sieve step by step and look at the dots: that is how many times a number was crossed out.

Up to a hundred the sieve makes 104 crossings, while there are 74 composites: thirty crossings are wasted. Switch to “Linear”, a different sieve described below: it crosses out each composite exactly once, by its smallest prime divisor, and makes exactly 74 writes. The label under a number shows how the linear sieve decomposed it: 60 = 2 · 30.

Why exactly the primes remain

The same idea gives the stop at √n: a composite number up to n always has a prime divisor at most √n. Trial division of a single number stops at the square root for the same reason.

Code

An array of “is it prime” flags, a byte per number. In Python crossing out the multiples of one prime is a single slice assignment, is_prime[p*p::p] = bytes(count): the loop runs inside the interpreter, in C, and that is tens of times faster than a Python for. In Go we write an ordinary loop: the compiler makes the same of it. The program prints the primes up to a hundred and counts how many crossings the sieve makes for various n.

         n   primes  crossings  per number  ln ln n
      1000      168       1411       1.411    1.933
   1000000    78498    2122048       2.122    2.626
  10000000   664579   22850051       2.285    2.780

Look at the “per number” column: from a thousand to ten million n grew ten thousand times, while the average number of crossings went from 1.4 to 2.3. The sieve is almost linear, and the next section explains where the “almost” comes from.

What it costs

A prime p crosses out about n/p numbers. The total number of crossings is less than

The sum of the reciprocals of the primes grows as a double logarithm: this is Mertens’s theorem of 1874, and the constant 0.2615 is his too. The double logarithm grows so slowly that it is hard to tell from a constant: for n = 10⁷ it equals 2.78, for n = 10¹⁸ 3.72. Hence O(n log log n) time. Memory is O(n): a byte per number, or, if you store bits and only odd numbers, sixteen numbers per byte.

Trial division of each number up to n is much more expensive: a prime k needs about √k divisions, giving on the order of n√n/ln n in all. Below is a measurement: four methods on the same n, with logarithmic axes, so the slope of a line shows the power. Trial division has a slope of about 1.4, one and a half powers slightly pressed down by the logarithm; the sieves about 1.

Two things on the chart are not obvious. The linear sieve, with the fewest writes, is no faster than the ordinary one: in Go they run neck and neck, and in Python the linear one lags about forty times, since its inner loop is written in Python, while the ordinary sieve’s loop is hidden in a slice. And at the largest n the plain sieve in Go gets more expensive faster than the segmented one. That is not the algorithm but memory. A 100 MB array fits in no processor cache, and a step with a large p lands in a new line of memory every time. The segmented sieve does the same crossings in pieces that fit in the cache, and at 10⁸ wins almost threefold.

Now the same in your browser. Choose n and run: the methods go one after another, and the table shows memory, number of operations and time.

At a billion only the segmented sieve remains: a byte per number is a gigabyte, while the segmented one needs a 32 KB window and a little over three thousand primes up to √10⁹ ≈ 31,623. The answer 50,847,534 agrees with the tables.

The linear sieve

Can each composite be crossed out exactly once? Yes, if it is assigned to one particular prime, its smallest divisor. Each composite m is written uniquely as m = p · i, where p is the smallest prime divisor of m; then all the prime divisors of i are at least p. So it is enough, for each i, to go through the primes p no greater than the smallest divisor of i and record p as the smallest divisor of the number p · i. Each pair (p, i) gives its own number, and each composite comes from exactly one pair: there are as many writes as composites. This is the linear sieve, whose O(n) running time was shown by David Gries and Jayadev Misra in 1978.

The program confirms it: up to 400,000 there are 366,139 composites, and exactly as many writes. But the main value of the linear sieve is not speed (above it is clear that in practice it does not give any) but a by-product: the table of smallest prime divisors spf. With it any number up to n is factored in O(log n) divisions: divide by spf[m] until one remains. That is how the program found the number of divisors of every number up to 400,000: 332,640 has the most, a whole 192. The price is memory: four bytes per number instead of one bit.

A table of smallest divisors can also be obtained with the ordinary sieve: when crossing out, write p into the number if it is still empty there. The repeated crossings remain, but the loop is simpler, and in practice the difference is small.

The segmented sieve: a window anywhere

To cross out the composites in a window [L, R), no array from zero is needed: the primes up to √R and an array the size of the window itself are enough. For each prime p find the first multiple in the window, ⌈L/p⌉ · p, but no less than p² so as not to cross out p itself, and step through the window with step p. This is the segmented sieve. Cutting the range from 0 to n into 32 KB windows gives the same sieve with O(√n) memory, and each window lives entirely in the cache, hence the gain in the measurement above.

The window can lie far from zero. The program below counts the primes up to ten million in windows of 32,768 numbers, and then sieves a million numbers after 10¹²: for that the 78,498 primes up to a million are enough, while the whole array from zero would take a terabyte.

[10^12, 10^12 + 10^6): 36249 primes, using 78498 primes up to 10^6
first 1000000000039, last 1000000999999
widest gap 236 after 1000000494263

There are 36,249 primes in this million, about one per 27.6 numbers, just as the prime number theorem promises: ln 10¹² ≈ 27.6. And the widest stretch without primes is 236 numbers in a row.

The article’s first scene is a segmented sieve too, only an incomplete one. For a 309-digit x no sieve reaches √x: that is a 155-digit number. So the window is sieved only with small primes, up to a thousand or up to a million. The survivors are not necessarily prime; they are numbers without small divisors, and the Miller–Rabin test does the rest. The share of survivors is predictable too, by Mertens’s third theorem:

where γ ≈ 0.5772 is Euler’s constant. For B = 719 the formula gives 17.1%, the exact product 17.0%, and in the window from the listing below 344 candidates of 2000 survived, 17.2%. The denominator is a logarithm, so each next tenfold growth of B cuts less and less: from a thousand to a million the share falls only by half, from 16% to 8%, while the divisions of the huge number by small primes become 470 times more.

The number here is “random” only in pretend: SHA-256 blocks of a phrase, so that Python and Go find the same one. From the start to the prime there are 457 odd candidates. A sieve with the 127 primes up to 719 cut the expensive checks to 76, a sieve up to a hundred thousand to 44.

Where it works

• Keys. Before the Miller–Rabin test OpenSSL weeds out candidates with the first primes: 127 odd ones up to 719 for 1024-bit numbers and 383 up to 2657 for 2048-bit ones. In the function BN_generate_prime_ex the remainders of the random start on division by these primes are computed once, then the candidate is shifted by 2, and the shifts at which a remainder becomes zero are discarded: the same sieve on a window, only lazy.
• primesieve. Kim Walisch’s open-source library for fast prime generation: a segmented sieve of Eratosthenes with wheel factorisation (the multiples of the very smallest primes are dropped in advance) and a window size tuned to the L1 and L2 caches of the particular processor.
• Goldbach’s conjecture. Tomás Oliveira e Silva, Siegfried Herzog and Silvio Pardi checked that every even number up to 4·10¹⁸ is a sum of two primes. The primes for this came from a cache-friendly segmented sieve; it took about 782 CPU-years. The work was published in 2014.
• Competitive programming. A table of smallest divisors up to 10⁶–10⁷ is a standard preparation: factoring numbers, Euler’s function, the number of divisors for many queries at once.

Where the sieve is not needed and where people trip on it

• A single number. To check one number for primality, no sieve is needed. For a 12-digit one, trial division up to √n, a million divisions, is enough, and for 20-digit and larger ones the Miller–Rabin test, a few exponentiations. The sieve pays off when all the primes in a row are needed, or a check of many nearby numbers.
• Memory. A byte per number up to 10⁹ is a gigabyte. Bits and odd numbers only give 62.5 MB, segments tens of kilobytes.
• Overflow. p * p does not fit in a 32-bit integer for p ≥ 46,341. In Go int is 64-bit and Python integers do not overflow, but in C and Java people trip on this constantly.
• A window with one. If a window starts at 0 or 1, they must be crossed out separately, and primes in the window itself must not be crossed out: hence “no less than p²” in the segmented sieve.

History

The writings of Eratosthenes, who lived in the 3rd century BC, have hardly survived. We know of the sieve from the “Introduction to Arithmetic” of Nicomachus, written around AD 100: Nicomachus attributes the method to Eratosthenes and describes it on odd numbers only, so the “odd only” optimisation is older than the word “algorithm” itself. The linear sieve with a proof of linear time was published by David Gries and Jayadev Misra in 1978. The segmented sieve has been known since the 1970s; one of the first works is a 1977 paper by Carter Bays and Richard Hudson, who sieved primes up to 10¹². In 1981 Paul Pritchard built a wheel sieve that needs even fewer than O(n) operations, but in practice a well-written sieve of Eratosthenes with a wheel is no worse.

Problems

The first is the segmented sieve in pure form, on windows far from zero. In the second the sieve is needed not for primes: first you will have to work out what to count.