Visualizing What is and What is Not the Sieve of Eratosthenes

The Classic Algorithm

The Sieve is a classic algorithm to create a list of prime numbers. However, this algorithm is often presented wrong, as discussed in an absolutely delightful paper by Melissa O'Neill. This is especially true when people present functional implementations of the algorithm. This post will explore how we can visualize the differences between implementations.

So what is the classic algorithm? To find all the prime numbers less than or equal to a given integer n, Wikipedia gives us the following method:

1. Create a list of integers from 2 through n: (2, 3, 4, ..., n).
2. Let p equal 2, the first prime number.
3. Enumerate all multiples of p by counting to n in increments of p, and mark them as composite.
4. Let p' equal p and let p equal the first number greater than p' in the list that is not marked composite. If there is no such number, stop. Otherwise, repeat the algorithm from step 3.

Okay, great, what does this look in actual code then? Well, I'm going to implement it using Java and Processing. I'm doing this because I think it'll be way easier to visualize the algorithm this way.

void classic_seive(int n){
  boolean[] l = new boolean[n + 1];
  // l is a table denoting whether an integer is composite 
  int p = 2;
  while(p <= n){
    for(int marker = 2 * p; marker <= n; marker += p){
      assert(marker % p == 0);
      l[marker] = true;
      visit(marker);
    }
    for(p++; p <= n && l[p]; p++){
      // increment p until l[p] is false.
    }
  }
}

This is a pretty straight-forward implementation of the algorithm. The only real note-worthy addition is the visit() call. All this does is tell my visualizer to register that value in the visualization. So what does this algorithm look like? I think the best way to capture this sieve is to look at the numbers being checked over time and this is exactly what the visit() call is going to do for us.

That graph shows which numbers are being checked over time as the sieve computes all primes less than 1000. If you look at what's happening in that graph, you'll see that early on we're eliminating lots of composite numbers and that later primes are eliminating fewer and fewer (this makes sense, as you only check the multiples of the prime and the gaps between multiples keep getting larger.

The Often Presented Functional Alternative

One problem with the sieve is that the algorithm is incredibly imperative. You initialize a large list and then mutate it a whole bunch as you repeatedly iterate through it. This does not lend itself to clever or beautiful functional implementations. However, people often supply an example implementation which is really quite lovely. From SICP:

(define (sieve stream)
  (cons-stream
   (stream-car stream)
   (sieve (stream-filter
           (lambda (x)
             (not (divisible? x (stream-car stream))))
           (stream-cdr stream)))))

I love the way Scheme reads. Unfortunately, this is not a true implementation of the sieve! Examine the recursive operation. The filter is applied to all remaining elements in the stream. Of course lazy evaluation muddles this a bit, but the whole cdr is being filtered. Every number still in the sieve and which is greater than the current prime will be checked for divisibility. Therefore, every composite number m will be checked once for every prime less than it's first prime factor. In the true sieve, every composite is checked once for every prime factor of that composite number.

Let's visualize this difference by reimplementing (a non-lazy version) in Processing:

void func_seive(int n){
  ArrayList vals = new ArrayList();
  for (int t = 2; t <= n; t++){
    vals.add(t);
  }
  vals = func_seive_helper(vals);
}

ArrayList func_seive_helper(ArrayList vals){
  if(vals.size() == 0)
    return vals;
  int cur = vals.get(0);
  ArrayList filtered = new ArrayList();
  for(Integer i : vals){
    if(i == cur){
      continue;
    }else{
      visit(i);
      if(i % cur != 0){
        filtered.add(i);
      }
    }
  }
  filtered = func_seive_helper(filtered);
  filtered.add(0, cur);
  return filtered;
}

My implementation is significantly uglier. I blame this on my comparably bad coding techniques (I cannot compete with SICP) and my need to reimplement filter. To further illustrate the difference, I changed visit so that the red hue of the points change the more that a particular number is checked. I actually did this before, but because numbers are only checked at most 3 times in the true sieve, the hue shift is basically imperceptible.

As you can see, numbers are being visited way more often. For example, 11 is visited when p equals 2, 3, 5, and 7. In the true sieve, 11 is never visited!

But What About Incremental Construction?

The functional "sieve" is not without its merits. In the classic sieve algorithm, you have to commit to the size of the sieve at the start of the algorithm. In the functional sieve, if you use lazy computation or streams then you can construct arbitrarily sized sieves. That's a nice difference, especially if you don't know at the outset how many primes you need to generate.

Of course, you can achieve a similar property by some hacking on the imperative classic sieve. Instead of constructing a sieve at whatever size you want initially, you use some "starter" size, and as the sieve needs to get larger, you resize your table, and update any new entries by rechecking them with the previously computed primes. In order to do this, your table needs to be a bit more complicated. This is because you need to remember what the last multiple was for each prime. That is easy enough to do by just upgrading our boolean table to an int table.

void fix_table(int cur_max, int[] l){
  int p = 2;
  while (p <= cur_max){
    if(l[p] == 0){
      l[p] = p;
    }
    for(int marker = p + l[p]; marker <= cur_max; marker += p){
      assert(marker % p == 0);
      l[marker] = 1;
      visit(marker);
      l[p] = marker;
    }
    for(p++; p <= cur_max && l[p] == 1; p++){
      // increment p until l[p] is false.
    }
  }
}

void incremental_seive(int n){
  int[] l = new int[n + 1];
  int table_size = 100;
  int p = 2;
  fix_table(table_size, l);
  while(table_size < n){
    int prev_size = table_size;
    table_size *= 2;
    if(table_size > n){
      table_size = n;
    }
    l = Arrays.copyOf(l, table_size + 1);
    fix_table(table_size, l);
  }  
}

So what does this one look like?

Pretty cool looking, I think. Also, notice that no numbers are "turning red" in this implementation, either.