The Sieve of Eratosthenes

The Ancient Greeks never fail to amaze me. Sometimes their ideas border the merely superstitious, sometimes you find a guy like Eratosthenes (c. 275 BC – c. 195 BC) that comes up with a truly modern idea. I’m not talking about his invention of geography (in the modern sense of the word) nor about his suspiciously accurate estimate of Earth’s size. No, I’m talking about an algorithm that has a very modern feel to it: the sieve of Eratosthenes

320px-Sieve_(PSF)

So what is it that I find so modern about this algorithm? Well, it has a computational feel to it, replacing smarts with a simple procedure, regular and repetitive. So, how hard is Eratosthenes’ sieve to implement? Not every: about 10 lines of code is needed in C++ (and possibly quite fewer if we use some other language).

The basic sieve of Eratosthenes is as follows. We start with the natural numbers laid out like this.

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 …

The first prime number is 2, so we cross every other number, all multiples of 2:

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 …

The next prime number, the next non-crossed number, is 3. We do the same for the multiples of 3:

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 …

Then we do the same for the next non-crossed number, 5, and so on. Continuing this procedure for ever, we find all prime numbers. Of course, here we have a rather short list, and so will you when you try it out, but you have to imagine that ‘…’ means forever.

*
* *

Neglecting the difficulty that integers are infinitely many and unbounded in magnitude, a 32bits type-safe implementation of the sieve could be as follows:


int main()
 {
  // scary, but only the shell
  // of the bitmap is on the stack
  std::vector<bool>
   erastobits(std::numeric_limits<uint32_t>::max());

  const uint32_t sqrt_max 
   = std::numeric_limits<uint32_t>::max() >> (8*(sizeof(uint32_t)/2));

  erastobits[2]=true;
  uint32_t lb=2,b=3;

  while (b<=sqrt_max)
   {
    std::cout << b << std::endl;

    // mark
    for (uint32_t l=0,c=b;l<c;l=c,c+=b)
     erastobits[c]=true;
    
    // find next
    while ( (lb<b) && erastobits[b]) lb=b,b+=2;
   }

  // printout the rest
  while (lb<b)
   {
    while ( (lb<b) && erastobits[b]) lb=b,b+=2;
    if (lb<b)
     {
      std::cout << b << std::endl;
      lb=b;b+=2;
     }
   }

  return 0;
 }

At first, we might remark that creating a vector of 2^{32} booleans on the stack is reckless, but 1) only the stub of the object is on the stack, the array is allocated on the heap and 2) std::vector<bool> has an explicit specialization that allows to use exactly one bit per boolean, thus reducing the needed memory from 4GB to merely 512MB (which still seemed quite a luxury not so long ago).

This program enumerates all 203280221 primes under 2^{32} in roughly one minutes on a 3.something-GHz i7 processor. Enumeration by factorization would take much longer.

Note the curious for condition lb<b. That’s one of the ways I have found to deal with a complete scan of (unsigned) integers. What happens is that eventually b+=2 overflows so that the result is (b+2) % 4294967296 (if you’re using 32 bits), and since 2 is always small enough compared to 4294967296; the for condition detects that the counter wrapped around, and the for terminates.

Note also that I test only until \sqrt{2^{32}} because all numbers beyond this point have either a larger divisor (in which case we already found the smallest) or a small divisor (which we already marked out). Therefore, testing until \sqrt{2^{32}} is enough for the sieve; and the last part of the program just printing out the last primes.

*
* *

It always come as a surprise to find some truly sophisticated piece of thinking coming from more than 2000 years ago. Yet, Eratosthenes lived in the same Era as Euclid, and a bit before the Antikythera device, so by now we should in fact ask ourselves just how much this civilization actually knew.

One Response to The Sieve of Eratosthenes

  1. […] we generate the prime numbers from 2 to . This could be done, for example, as we did before. Once the list is acquired, we will proceed to test each of those prime numbers as the , and see […]

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: