Fast Prime Calculation - Part 1 : Sieve of Eratosthenes

View / download this post as a PDF

When working on Project Euler problems, you quickly reach a point where you want a fast way to compute primes. Note: I prefer to have my solutions compute a list of primes rather than loading a list of precomputed primes from a file, since I use my solutions as a test bed for my numerics library. In this post, I'll outline a method similar to the Sieve of Eratosthenes and then in a later post I'll go over another iteration of this method that can work within a limited memory budget as well as giving it a significant speed up.

Basic Prime Computation

First off, what is a prime number? A prime number is a natural number (positive, non-zero number) that has exactly 2 divisors: 1 and itself (note: 1 is not considered a prime number). The simplest test you might think of to see if a number n is prime is to divide it by all smaller numbers to make sure none of them divide evenly into n, and if so, n is prime.

How can we make this faster? Let's prove that you only need to check potential divisors up to √n. By contradiction, assume that you might find your first divisor d₀ > √n. If that's the case then you get another divisor by dividing n / d₀ = d₁. Either d₁ is greater than or equal to d₀ in which case d₀ · d₁ > √n · √n and so d₀ · d₁ > n, which is a contradiction. Otherwise d₁ < d₀ in which case we should have discovered it earlier, which contradicts our assumption that d₀ is the first divisor we find. Either case is a contradiction, so we know that if no divisors are found less than or equal to √n, then n is prime.

Now, let's assume you find a number d₀ that divides n. Either d₀ itself is prime, or there is another smaller number d₁ that divides d₀ evenly. Continuing in this fashion, you can break n down into divisors, which can be broken down into their own divisors, on and on until you're left with a product of prime numbers. This is known as the prime factorization of n. Given this fact, you don't have to check all numbers up to √n, instead you just need to check all the primes less than or equal to √n. From this we can see that if you want to generate a list of prime numbers, it's logical to start with the smallest primes (2, 3, 5, ...) and gradually work your way up, using your previously computed primes to check the new candidate numbers.

Skipping Multiples

Sometimes the best way to see how you can speed a process up is to do it by hand. Let's take a look at what it's like to test the first few numbers.



We quickly see that 2 and 3 are prime. Then we see that 2 divides 4, so 4 isn't prime. Then, since 2 doesn't evenly divide 5 (3 is greater than √5 so we don't have to check it), we see that 5 is prime. Again we see that 2 divides 6, ruling it out. Then 7 is prime. Again, 2 divides 8, ruling it out. We start to see that we can easily skip all the even numbers, since we know they will always have 2 as a divisor.

What if we want to do the same thing for multiples of 3. If we look at all the numbers that are multiples of 2 and/or 3, they form a pattern that repeats every 6 numbers.



We can see that in every group of 6 numbers, the first and the fifth are the only numbers we need to consider for our prime list, since all the other numbers will either have 2 or 3 as a divisor. In fact this process can be extended to any number of primes that you want to rule out. Just take a list of primes p₀, p₁, p₂, ... and make a table with size p₀ · p₁ · p₂ · ... and then mark off all the multiples of each prime. The remaining numbers are your potential primes, and the pattern of potential primes can be repeated to rule out further non-prime numbers.

Quick Iteration Through Candidates

Again, if we imagine doing this process manually, you still need to search through the table to find the next potential prime slot. We can optimize this a bit by turning the table into a list of offsets that allow us to easily and instantly step over consecutive non-primes to take us from one candidate to the next. The following figure shows the table generated with 2, 3 and 5 (size 30) and the corresponding offsets.



Our array of offsets is (6, 4, 2, 4, 2, 4, 6, 2). Note that we cross out 2, 3 and 5 since, upon repeating the pattern, those slots correspond to numbers that are multiples of those primes and so they need to be skipped. Also note the last candidate is 29, and in order to loop back around to the beginning of the table, the last offset is 2. It turns out this last offset to restart the pattern is always 2. Let's take a moment to confirm this fact. The last number in a table is obviously divisible by all the primes used to generate the table, since it is, by definition, the product of all of those primes. Thus, one number less than that can not be divisible by any of those primes, since the smallest possible prime is 2 so you would need to walk at least 2 slots away from the last slot before any prime could divide it evenly. Similarly, the first slot will always be a candidate since it's only 1 slot away from the end of the table (upon repeating the pattern). Therefore since slot size - 1 is always a candidate and slot 1 is always a candidate, the last offset is 2 in order to loop back around to the beginning of the table.

Building Offset Lists

Now let's consider how to build this list of offsets. We could build the table and then walk over it, subtracting positions of candidate slots, but let's consider a method that doesn't need the table at all. If we start with the offset list for a table generated with only the prime number 2, the list has only one offset with a value of 2. This matches the conclusion we made in the previous paragraph that the last offset is always 2, as well as matching our intuition that we need to keep adding 2 to skip over all the even numbers.

Now if we want to grow this list to account for the multiples of 3 as well, we should repeat the list 3 times (so that it covers a table of size 6). Then we can walk through the offsets adding them up and any time we land on a multiple of 3 we can combine that offset with the next one to make us skip that candidate. So the offset list for the table generated by 2 and 3 is (4, 2), as shown in the following figure



This can be repeated for as many primes as you like.

Once you grow your offset list to an acceptable size, you can use it to walk through candidate primes as far as you want, and for each candidate, check it against all primes up to its square root (remembering that you don't need to try any of the primes used to generate your offset list, since you know they won't divide evenly).

Offset List Size

Here we'll show by induction that the number of offsets k_n in a list generated with primes p₀, p₁, ... , p_n is



Before we start, let's establish one preliminary fact. When you're dealing with the integers modulo a prime number p, you can take any number a from 1 to p-1 and list its multiples (modulo p) out to p · a, and you'll get all the numbers from 0 to p - 1 with no repetitions. This is due to the fact that the set of integers modulo a prime is a simple cyclic group under addition. In particular, note that the last value generated is 0, because p · a is divisible by p. In other words, if you list all the multiples up to (p - 1) · a then you'll get all the numbers from 1 to p - 1 (the same as the above statement with 0 excluded).

Getting back to our table and list of offsets, if we only have one prime p₀, then the table size is p₀ and the only number marked out will be p₀ itself. So it would have p₀ - 2 offsets (each with a value of 1) before p₀ and then the final offset of 2 to loop back to the beginning giving a total of p₀ - 1 offsets. This establishes the base case for our inductive proof.

Next assume that for a list of primes p₀, p₁, ... , p_n-1 the number of offsets generated is



First note that the table size for the set of primes up to p_n-1 is



Now let's introduce the next prime p_n. We repeat the offset list from the previous step p_n times and then proceed to combine offsets that would have taken us to a spot in the table that is a multiple of p_n. Now consider any of the single offsets from the previous list. Let's call it f_i, and let's call the slot in the table to which this offset leads s_i. Also we'll say that s_i = x_i mod p_n. If x_i is zero, then that means that p_n divides s_i (and so the offset f_i should be combined with f_i+1). If x_i is non-zero, then it remains untouched

The previous offset list spanned a table of size t_n-1, and we should note that t_n-1 is not divisible by p_n since all of its factors are prime and none are equal to p_n (which is also prime). Therefore t_n-1 = m mod p_n and m is non-zero. From the previously stated fact about each non-zero member of the integers modulo p_n generating the entire group of integers modulo p_n, we know upon repeating the offset list p_n times, the duplicates of the offset f_i will lead to x_i, x_i + m, x_i + 2 m, ... x_i + p_n · m. In general these will be of the form x_i + j · m where j goes from 0 to p_n. For one of these, j · m will be equal to the additive inverse of x_i mod p_n and so x_i + j · m = 0 mod p_n.

Therefore every offset should have one and only one duplicate that should be cancelled out. Since there are k_n-1 offsets being duplicated p_n times and each being cancelled once, we can see that



which completes the proof.

One thing to notice from this fact is that if you have an offset list generated from a set of primes that doesn't include 2, repeating the process to include 2 will leave the number of offsets completely unchanged since for primes p₀, p₁, ... , p_n, the number of offsets is



Then adding p_n+1 = 2 to the list of primes you get



Again, this is because the offset list generated from the previous primes will be doubled in length, but then each offset will either be cancelled out in the first or second repetition. Therefore, it's always a good idea to include 2 in your list of primes, since it has no memory overhead. You can use similar logic to see that the smaller your primes are, the more effecient they are in terms of memory.

Conclusion

This method is very close to something known as the Sieve of Eratosthenes, where instead of deciding on a fixed number of primes, you start with an upper bound for the primes you want to compute, list all the numbers up to that upper bound and then you mark off multiples of each prime in turn. Once you've finished marking off all the multiples of a particular prime, the next prime is the first unmarked spot. This has the advantage of not needing to divide by a bunch of known primes to test candidates, since you immediately know the next prime after finishing marking off multiples of the previous prime, but it inherently requires you to put an upper bound on the primes you want to compute. Next time we'll look at another approach that's a little closer to the Sieve of Eratosthenes without needing the upper bound limitation and using some of what we've discussed in this post to help speed things up.