Fast Prime Calculation, Part 1: The Sieve of Eratosthenes Philip Nowell Co-founder, Controlled Panic LLC April 10, 2014 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.

Contents 1 Basic Prime Computation

2

2 Skipping Multiples

2

3 Quick Iteration Through Candidates

4

4 Building Offset Lists

5

5 Offset List Size

6

6 Conclusion

8

1

1

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 d0 > n. If that’s the case then you get another divisor by dividing√n/d√ 0 = d1 . Either d1 is greater than or equal to d0 in which case d0 ·d1 > n· n and so d0 ·d1 > n, which is a contradiction. Otherwise d1 < d0 in which case we should have discovered it earlier, which contradicts our assumption that d0 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 d0 that divides n. Either d0 itself is prime, or there is another smaller number d1 that divides d0 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.

2

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. 2

1 2 3 4 5 6 7 8 2

Factors

2

2

3

Figure 1: The first few numbers and the primes that evenly divide them 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. 1 2 3 4 5 6 2

2 3

7 8 9 10 11 12 2

3

2

2 3

13 14 15 16 17 18 2

3

2

2 3

19 20 21 22 23 24 2

3

2

2 3

Figure 2: Ruling out multiples of 2 and 3 From Figure 2, 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 p0 , p1 , p2 , ... and make a table with size p0 · p1 · p2 · ... and then mark off all the multiples of each prime. The remaining numbers are 3

your potential primes, and the pattern of potential primes can be repeated to rule out further non-prime numbers.

3

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. Figure 3 shows the table generated with 2, 3 and 5 (size 30) and the corresponding offsets. 6

4

2

4

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 2

4

6

2

16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 1

Figure 3: Ruling out multiples of 2, 3 and 5 and the offsets used to quickly step through the remaining numbers 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

4

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.

4

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) (see Figure 4)

2

2 5

3

1

2 1

4 2 1

2

2 5

3

1

Figure 4: Duplicating the offset list for a single prime (2) and then combining offsets when a multiple of 3 is encountered 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, 5

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).

5

Offset List Size

Here we’ll show by induction that the number of offsets kn in a list generated with primes p0 , p1 , ..., pn is kn =

n Y

pi − 1

i=0

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 p0 , then the table size is p0 and the only number marked out will be p0 itself. So it would have p0 − 2 offsets (each with a value of 1) before p0 and then the final offset of 2 to loop back to the beginning giving a total of p0 − 1 offsets. This establishes the base case for our inductive proof. Next assume that for a list of primes p0 , p1 , ..., pn−1 the number of offsets generated is kn−1 =

n−1 Y

pi − 1

i=0

First note that the table size for the set of primes up to pn−1 is tn−1 =

n−1 Y i=0

6

pi

Now let’s introduce the next prime pn . We repeat the offset list from the previous step pn times and then proceed to combine offsets that would have taken us to a spot in the table that is a multiple of pn . Now consider any of the single offsets from the previous list. Let’s call it fi , and let’s call the slot in the table to which this offset leads si . Also we’ll say that si ≡ xi mod pn . If xi is zero, then that means that pn divides si (and so the offset fi should be combined with fi+1 ). If xi is non-zero, then it remains untouched The previous offset list spanned a table of size tn−1 , and we should note that tn−1 is not divisible by pn since all of its factors are prime and none are equal to pn (which is also prime). Therefore tn−1 ≡ m mod pn and m is non-zero. From the previously stated fact about each non-zero member of the integers modulo pn generating the entire group of integers modulo pn , we know upon repeating the offset list pn times, the duplicates of the offset fi will lead to xi , xi + m, xi + 2m, ... xi + pn · m. In general these will be of the form xi + j · m where j goes from 0 to pn . For one of these, j · m will be equal to the additive inverse of xi mod pn and so xi + j · m ≡ 0 mod pn . Therefore every offset should have one and only one duplicate that should be cancelled out. Since there are kn−1 offsets being duplicated pn times and each being cancelled once, we can see that kn = kn−1 · (pn − 1) n Y = pi − 1 i=0

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 p0 , p1 , ..., pn , the number of offsets is kn =

n Y i=0

7

pi − 1

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

n+1 Y

pi − 1

i=0 n Y

! pi − 1

· (pn+1 − 1)

i=0

=

n Y

! pi − 1

· (2 − 1)

i=0

=

n Y

pi − 1

i=0

= kn 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.

6

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.

8