Faster Sieve of Eratosthenes
About a decade ago I was participating in a lot of programming challenges and one of those challenges was posted on a programming forum where people were interested in finding out who could make the most optimized version of common mathematical problems.
One of those common math problems is the ability to find out prime numbers. I took the challenge and implemented an optimized sieve of eratosthenes by assuming that the number 2 was a prime and halving the memory space required, like many others. My solution was the second fastest implementation of all the participants, according to the benchmarks, but both the first and second fastest solutions were essentially the same. But while the participant’s solutions were being benchmarked, I found a very discrete pattern that allowed me to come up with another solution that was, by far, faster than any solution presented in the challenge. I still got second place because the challenge “owner” (and winner of the challenge) didn’t want to benchmark my new solution.
Disclaimer: I’m not a mathematician. I’m a coder who happened to trip on this while participating on a challenge and that kept investigating up to a certain point.
My solution starts with the assumption that everyone knows that the numbers 2 and 3 are prime numbers. Everyone knows that, past the number 2, only odd numbers can be prime. But I took it one step further by eliminating all odd numbers that are multiple of 3 as well.
A list of points:
- Multiples of 3 are either odd or even (thanks Captain Obvious);
- Even multiples of 3 are composed by multiplying 2 and 3, so they’re essentially multiples of 6;
- There are no other odd numbers besides those before and after 6 that are not multiples of 3(because around odd multiples of 3 you’ll only find even numbers and with the odd numbers found around multiples of 6, you’re basically covering every number possible);
By doing this, you effectively cut a nice chunk of the memory requirements… except you can’t actually know how many memory slots to jump like this. Or can you?
The sieve is split into two columns: one that is 6*n-1 and other that is 6*n+1. Lets remove all multiples of 5 from the columns:
Do you see a pattern? It can be expressed as 6*(5*n+1)-1 and 6*(5*n–1)+1. Lets try to do the same for the number 7:
Notice the formulas become inverted. For this case we have to express is as 6*(7*n+1)+1 and 6*(7*n–1)-1. And the pattern gets messier when we move towards the second row of prime numbers:
This breaks our distance between non primes for the same prime across columns. Now the formulas for 11 are 6*(11*n+2)-1 and 6*(11*n–2)+1 and for 13 are 6*(13*n+2)+1 and 6*(13*n–2)-1.
But there’s still a pattern there. If you keep doing this for the next set of numbers, you’ll notice that the value that changes is incremented every time we move one row. E.g.: the formulas for 17 and 19 are 6*(17*n+3)-1, 6*(17*n–3)+1 and 6*(19*n+3)+1, 6*(19*n–3)-1 respectively. So, the change between these formulas are dependent on another variable: the index of the number in the column.
Knowing this, implementing a sieve that works the same way as the existing sieves becomes easy. But we can still make these formulas simpler. Lets call pc1 and pc2 the variables that correspond to prime numbers in the columns 1 and 2 respectively, and i the index (non zero based) where we can find the prime number. From what we know from before, this is what we have:
6*(pc1*n+i)-1, 6*(pc1*n-i)*n+1
6*(pc2*n+i)+1, 6*(pc2*n-i)-1
The i value can actually be calculated by applying a simple division by 6 on the closest multiple of 6 to the prime number (because its the nth multiple of 6 which this system is based on), which means that, respectively for each column:
i = (pc1+1)/6
i = (pc2–1)/6
So:
6*(pc1*n+(pc1+1)/6)-1, 6*(pc1*n-(pc1+1)/6)+1
6*(pc2*n+(pc2-1)/6)+1, 6*(pc2*n-(pc2–1)/6)-1
And because we should always simplify our formulas:
6*pc1*n+pc1, 6*pc1*n-pc1
6*pc2*n+pc2, 6*pc2*n-pc2
So, notice that now we have the same exact formulas for both columns, which means that we should discard the pc1 and pc2 nomenclature with something more generic such as p:
6*p*n+p, 6*p*n-p
We still have two formulas, one for the column matching the prime number and other for the opposite (respectively) which is enough to simplify the sieve implementation, but mathematically speaking, both formulas can be condensed into a single formula that yields two results:
6*p*n±p
Whatever the prime number p, this formula will always provide two non-prime numbers (for each iteration of n) next to multiples of 6 where p is one of its factors.
Using this strategy it’s not only possible to reduce the memory needs by 2/3 against a common sieve, but also save a lot of cycles and memory reads and writes.
You can find my implementation of this sieve @ https://github.com/KTachyon/FastPrimeSieve
The repo contains a XCode project, but you just need to compile the main.m file to get it running. On my machine it calculated the sieve for the first billion numbers in just under 4 seconds.
