# Incremental Sieve of Eratosthenes

The Sieve of Eratosthenes is a well known algorithm for computing primes, but suffers from space requirements. The Incremental Sieve solves that problem.

The Sieve works by walking through each of the numbers in the range. But that means that you have to store the numbers. Even if you reduce this to a bit map, the space requirement is O(n).

Thinking as a software engineer, it is evident that the core of the Sieve’s implementation would be a nested loop.

``````for each prime
for each multiple
mark number as composite
end
end
``````

The inner loop depends on the outer loop. That is good for efficiency, but not for analysis. So lets rewrite as:

``````for each prime
for each number
if number is a multiple of prime
mark number as composite
end
end
end
``````

Now we can see that we could get the same results by inverting the loops.

``````for each number
for each prime
if number is a multiple of prime
mark number as composite
end
end
end
``````

So, now, we don’t have to store the numbers. Instead we only have to store the primes. And since we are storing primes, we should make the test store it.

``````for each number
for each prime
if number is a multiple of prime
go to next number
end
end
add number to list of primes
end
``````

The remaining question is - how are we going implement that `if` ? We could just divide. But division is slow, so lets think of another way.

What if we stored a list of multiples for each prime ? But that quickly becomes O(n^^2) - until we realize that we really only care about multiples that are “near” our current target number. In, fact we only care about a multiple if it is equal or larger than our target. And we really only care about the smallest such multiple. So, in fact, we only need one multiple per prime. And so, the Incremental Sieve.

In pseudo-code it wold look like:

``````for each number
for each prime
loop while current multiple < number
get next multiple of the prime
end
if number == current multiple
go to next number
end
end
add number to list of primes
end
``````

Or in python:

``````#!/usr/bin/python3

primes = [ [2, 2] ]

for x in range(3, 100000) :
isprime = True;
for pm in primes :
while pm < x :
pm += pm
if pm == x :
isprime = False;
break
if isprime :
primes.append([x, x]);

print( len(primes))
``````

This takes about 21 seconds to run on my ARM Chromebook. Not bad, but maybe we could do better.

One thing to notice is that all primes (above 2) will be odd. So, there is really no reason to iterate through the even numbers. But that means we don’t need to store the 2 either, since it will never do us any good.

``````#!/usr/bin/python3

primes = [ [3, 3] ]

for x in range(5, 100000, 2) :
isprime = True;
for pm in primes :
while pm < x :
pm += pm
if pm == x :
isprime = False;
break
if isprime :
primes.append([x, x]);

print( len(primes))
``````

And that runs in about 20 secs. Hmmm. Not as helpful as I had hoped, but still worth it.

One fact to note about multiplication is that if `a * b = c` then one of those number must be below (or equal) to the square root of c and the other must be above (or equal) to the square root. So, if c is not prime, then it must have at least one of its factors below its square root. Which means, that if we can’t find a factor there, we can stop looking. For our top value of 100,000, we only need to search primes less than about 317. That is a huge savings.

``````#!/usr/bin/python3

import math

primes = [ [3, 3] ]

for x in range(5, 100000, 2) :
limit = int(math.sqrt(x))
isprime = True;
for pm in primes :
if pm > limit :
break;
while pm < x :
pm += pm
if pm == x :
isprime = False;
break
if isprime :
primes.append([x, x]);

print( len(primes))
``````

And this runs in about 0.6 seconds. Increasing the loop to 1,000,000, the progam takes about 9 seconds to run.

Just out of curiosity, I also implemented the last version in c++ :

``````#include <list>
#include <iostream>
#include <cmath>

using vtype = long long;

struct pm {
vtype prime;
vtype multiple;

pm(vtype p, vtype m) :prime(p), multiple(m) {}
};

std::list<pm> primes = { {3, 3} };

int main() {
for (vtype x = 5; x < 1'000'000; x+=2) {
vtype limit = std::sqrt(x) + 1;
bool isprime = true;
for (auto &[p, m] : primes) {
if (p > limit) break;
while (m < x) {
m += p;
}
if (m == x) {
isprime = false;
break;
}
}

if (isprime) {
primes.emplace_back(x, x*x);
}
}

std::cout << primes.size() << "\n";
}
``````

Compiling and running we get :

``````\$ g++ -std=c++17 -O3 incsieve.cpp
\$ time ./a.out
78497

real    0m0.093s
user    0m0.084s
sys     0m0.010s
``````

Well that is mighty nice. For the record, it took about 30 secs to run with the loop set to 100,000,000.

I hope you had fun looking at this variation on an old “favorite”.