1
\$\begingroup\$

The Sieve of Eratosthenes is one of the oldest known algorithms and yet it is still the best method for sieving prime numbers in bulk in many situations. It is brutally simple, but to make it brutally fast requires adding a quite few complications.

A lot of information about these complications is availabe here on Code Review, on Stack Overflow and elsewhere on the Web, but yet it seems surprisingly difficult to locate that information when one needs it, as it is scattered far and wide and intermingled with incredible amounts of noise. Hence I propose this series of articles which takes a simple, naive textbook implementation and introduces interesting complications one by one, including all the whys and wherefores. That way you can pick and choose when you craft your own version of the sieve for a specific purpose, or you can post your own version of a certain complication for others to study and admire.

I'll be using C# because it is readable for many programmers with other language backgrounds as well, and things like the built-in array bounds checking are big boon when doing rapid prototyping. Moreover, C# can be used with the free LINQPad which allows code & run cycles with a minimum of fuss and which adds a neat Dump() method to every class under the hood so that you can do things like

remaining_numbers(sieve_up_to(100)).Dump()

and see your results as a nicely formatted table. I usually use LINQPad to study algorithms and to prototype algorithmic optimisations, even if the final implementation language will be C++ or Delphi (or even T-SQL, in the case of check digit algorithms).

Here is a transcribed version of the algorithm given in the Wikipedia article Sieve of Eratosthenes:

static bool[] sieve_up_to (uint n)
{
    Trace.Assert(n < int.MaxValue);    // in .NET, size_t is signed integer (!)

    var eliminated = new bool[1 + n];  // arrays are 0-based but we will index by number

    eliminated[0] = true;  // 0 is not a prime, but the loop is unable to eliminate it
    eliminated[1] = true;  // 1 is not a prime, but the loop is unable to eliminate it

    for (uint i = 2, sqrt_n = (uint)Math.Sqrt(n); i <= sqrt_n; ++i)
        if (!eliminated[i])
            for (uint j = i * i; j <= n; j += i)
                eliminated[j] = true;

    return eliminated;
}

The outer loop can terminate once i exceeds the square root of n because every composite c ≤ n must have at least one prime factor p ≤ sqrt(n) and thus will have been crossed off during the round for that p.

The inner loop can start at i * i for an analogous reason: all composites c < i * i must already have been crossed off because for a given prime p the number p * p is the smallest composite that does not contain prime factors smaller than p.

The logic of the sieve array is the opposite of that in the Wikipedia pseudo code; that is, it contains true for crossed-off numbers (non-primes).

There are three reasons for this. First, the old Greek didn't say 'make a cross next to all numbers that you want to sieve and then start removing the crosses in the following manner...'. Second, it is more logical. We're trying hard to something unto all composites and at the end we look at what's left unscathed; hence it makes sense to record what we're doing unto composites, instead of undoing a superfluous operation that was applied undiscriminately to all numbers in the range. Third, it is both less error-prone and more efficient to rely on the implicit 0/false/null initialisation guaranteed by the runtime, especially as it cannot be suppressed anyway. (Sidenote: For more advanced sieves it can make sense to use different conventions, in order facilitate efficient decoding.)

The nested loop structure shows that the algorithm operates on two levels: on the outer level it enumerates the sieving primes, and on the inner level it does the actual sieving which - among other things - also produces new sieving primes. This is possible because the sieving range is anchored at the bottom and thus also includes any sieving primes that may be needed.

If the job of obtaining the sieving primes is delegated to a different function then the algorithm can be written like this, to make its structure and logic more apparent:

foreach (var prime in primes_up_to((uint)Math.Sqrt(n)))
    for (uint i = prime * prime; i <= n; i += prime)
        eliminated[i] = true;

Of course, this offers no advantage if there is only ever one sieving to be done during a program run. It is, however, the first step towards sieving ranges that are not anchored at zero and which might therefore not include the necessary sieving primes.

The ability to sieve an arbitrary interval [m, n] instead of a zero-anchored range [0, n] can be helpful in several contexts - like when partitioning a range for the distribution of work in time or across several processors, or when solving coding challenges like the SPOJ problem PRIME1.

PRIME1 asks for the primes in a window that can go as high as 1,000,000,000 but which is at most 100,000 numbers wide. In such a case, sieving all primes up to the upper end of the window instead of just the primes in the target window would mean doing up to 10,000 times the amount of work that's necessary. A well-written sieve could still beat the generous time limit of 6 seconds but a windowed sieve takes only a fraction of a millisecond, a speed gain of four orders of magnitude.

Windowed sieving adds the complication of having to compute the point where the crossing-off sequence for a prime p first intersects with the sieving window. There are two cases:

1) If the starting point p * p of the striding sequence lies at or above the start of the window then this is as simple as subtracting the window base.

2) If the starting point lies before the start of the window then a bit of modulo magic is required. With window base W, starting point K and stride S, the expression S - (W - K) % S gives the right result except for the the special case (W - K) % S == 0. In that case the desired result is 0 but the expression gives S. This can be fixed by aiming the expression at W - 1 instead of W, since adding a value between 1 and S to W - 1 results in a value between W + 0 and W + S - 1 as desired. Hence the expression becomes

(S - 1) - (W - 1 - K) % S

Et voilà, we can now sieve an arbitrary range of numbers:

static bool[] sieve_between (uint m, uint n)
{
    Trace.Assert(n - m < int.MaxValue);

    uint sieve_bits = n - m + 1;
    var eliminated = new bool[sieve_bits];

    for (uint i = m; i < 2; ++i)
        eliminated[i] = true;

    foreach (uint prime in primes_up_to((uint)Math.Sqrt(n)))
    {
        uint start = prime * prime, stride = prime;

        if (start >= m)
            start -= m;
        else
            start = (stride - 1) - (m - start - 1) % stride;

        for (uint j = start; j < sieve_bits; j += stride)
            eliminated[j] = true;
    }

    return eliminated;
}

There is a separate variable stride for the step width because this value does not have to be equal to the current prime.

For example, for all primes p greater than 2 the number p * p + p must necessarily be even and hence must already have been crossed off. The same goes for all other slots whose distance from p * p is an odd multiple of p. This means we can skip the even multiples by setting

stride = prime == 2 ? prime : 2 * prime;

for all primes > 2 and so save quite a bit of work. One simple way of accomplishing this without a conditional is by using the following initialising expression:

stride = prime << (byte)(prime & 1);  // cast necessary because of a C# language quirk

This is the simplest form of 'wheeled striding' - i.e. skipping multiples of small primes when stepping through the sieve array. IOW, we step only on slots where the corresponding wheel has a potentially prime-bearing spoke. For the mod 2 wheel, only the odd spoke can bear primes whereas the even spoke is always composite, with the exception of the slot for the wheel prime 2 itself.

Using wheeled striding in a sieve that is not wheeled or has a lesser wheel order realises only some of the possible benefits. In the example here we still had to allocate space for all even numbers, and the even numbers still had to be crossed off during the first round of the outer loop (the one for the prime 2). But we did halve the number of all subsequent crossings-off with just a minor index calculation trick, and this resulted in a noticeable speed-up. Wheeled striding can also be used profitably when enumerating ranges of numbers in a loop, like in a divisibility pre-check for big numbers.

Here are the timings for sieving a window of the given width whose lower end is 0, without and with the double stride trick:

    window                            stride = prime       stride = p << (p & 1)
--------------------------------------------------------------------------------
     32768 @ 0:     3514 primes       0,0 ms 1043,6 k/ms      0,0 ms 1588,1 k/ms
    100000 @ 0:     9594 primes       0,2 ms  657,9 k/ms      0,1 ms  812,1 k/ms
   1000000 @ 0:    78500 primes       2,5 ms  407,1 k/ms      1,6 ms  643,3 k/ms
  10000000 @ 0:   664581 primes      36,1 ms  277,3 k/ms     25,0 ms  400,7 k/ms
 100000000 @ 0:  5761457 primes     835,8 ms  119,6 k/ms    518,9 ms  192,7 k/ms
1000000000 @ 0: 50847536 primes    9995,7 ms  100,0 k/ms   5947,3 ms  168,1 k/ms

The timings show clearly that the simple sleight-of-hand trick with stride results in a noticeable speed-up.

Another option for the double stride thing would be to split out the round for the prime 2 and let the for loop start with the prime 3, so that all rounds can use the same value for the stride. The round for the prime 2 could then be replaced with blasting the pattern true, false, true, false, ... over the sieve with some buffer copy magic and adjusting for the prime 2 itself afterwards.

This is a simple example of the technique called 'presieving', whose speed gain comes from avoiding the crossing-off work for the handful for tiny primes that have the most 'hopping' to do because their strides ares so tiny. However, it is absolutely not worth it doing this here - there will be a separate article for this after we have harvested all the lower-hanging fruit. Besides, the sleight-of-hand trick with even multiples will become superfluous when odds-only sieving is introduced in the next article of this series.

And now the timings for sieving windows whose upper end is 10^9. This means that all sievings are using the exactly same number of sieve primes, and this in turn means that the speed differences are mostly due to the relation between the sieve size and the size of the memory caches of the CPU:

    window                                               stride = prime       stride = p << (p & 1)
---------------------------------------------------------------------------------------------------
     32768 just before 1000000000:     1572 primes       0,1 ms  481,0 k/ms      0,0 ms  667,9 k/ms
    100000 just before 1000000000:     4832 primes       0,3 ms  386,8 k/ms      0,2 ms  588,1 k/ms
   1000000 just before 1000000000:    47957 primes       3,4 ms  296,0 k/ms      2,0 ms  495,8 k/ms
  10000000 just before 1000000000:   482825 primes      57,2 ms  174,9 k/ms     40,8 ms  245,4 k/ms
 100000000 just before 1000000000:  4838319 primes    1023,5 ms   97,7 k/ms    621,1 ms  161,0 k/ms
1000000000 just before 1000000000: 50847536 primes    9965,5 ms  100,3 k/ms   6080,9 ms  164,4 k/ms

As you can see, the sieving speed decreases as the sieve size exceeds first the L1 cache (32 KiB) and then the L2 cache (256 KiB), and it finally plateaus after significantly exceeding the L3 cache size of 8 MiB. Optimising the sieve for efficient cache usage is one goal of segmented sieving, which will be dealt with in a later article of this series.

Please observe the second line, which deals with the worst case for the PRIME1 scenario. As promised, just fractions of a millisecond.

Sidenote: apart from some special cases it can be advantageous to leave the sieving part and the part that extracts the primes from the sieve separate. For one thing, the separation facilitates unit testing and experimentation with various optimised implementations. For another, the specific task at hand may not require the actual primes at all but rather their count or sum or digest, or it may ask only for every thousandth prime. So you can leave the actual sieving code unchanged and simply write a processing function tailored to the task at hand.

Here is a little helper function for reading out the numbers that haven't been crossed off:

static List<uint> remaining_numbers (bool[] eliminated, uint sieve_base = 0)
{
    var result = new List<uint>();

    for (uint i = 0, e = (uint)eliminated.Length; i < e; ++i)
        if (!eliminated[i])
            result.Add(sieve_base + i);

    return result;
}

This can be used to obtain the primes in a given range [m, n] without having to compute all the primes up to n:

static List<uint> primes_between (uint m, uint n)
{
    return remaining_numbers(sieve_between(m, n), m);
}

That's it, the foundations of sieving à la Eratosthenes are all there. All that's left to do for further instalments is to take this poor code and disfigure it by adding one complication after the other - just for the sake of making it faster by orders of magnitude...

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.