EPI: Optimized Algorithm for Generating Primes Numbers

Question

I am going through the "Elements of Programming Interview" in python currently and have gotten stuck at this part. The code below generates primes up to n. The explanation is rather lacking and i do not have the mathematical background to make sense of it.

We can improve run-time by sieving p's multiples from p^2 instead of p, since all numbers of the form kp, where k < p have already been sieved out.

The code is below:

def generate_primesII(n):

    if n < 2:
        return []

    size = (n - 3) // 2 + 1
    primes = [2]  # stores the primes from 1 to n

    # is_prime[i] represents (2i + 2) is prime or not
    # Initially set each to true. Then use sieving to eliminate nonprimes
    is_prime = [True] * size

    for i in range(size):
        if is_prime[i]:

            p = i * 2 + 3
            primes.append(p)

            # Sieving from p^2, where p^2 = (4i^2 + 12i + 9). The index in is_prime
            # is (2i^2 + 6i + 3) because is_prime[i] represents 2i + 3

            # note that we need to use long for j because p^2 might overflow
            for j in range(2 * i**2 + 6 * i + 3, size, p):
                is_prime[j] = False
    return primes

My questions are:

1. how did they come up with the formula for the size
2. they say is_prime[i] represents (2i + 3) is prime or not. I cant make sense of why 2i + 3.
3. how did they get p = i * 2 + 3
4. what does the following mean Sieving from p^2, where p^2 = (4i^2 + 12i + 9). The index in is_prime is (2i^2 + 6i + 3) because is_prime[i] represents 2i + 3
5. why does the range of j begin with 2 * i**2 + 6 * i + 3

Most of the numbers seem rather random to me

Answer 1

There are two key tricks which are simultaneously done here. That, I believe, is the main reason behind your confusion. The first is a mathematical fact about the progression about the sieve algorithm. (i.e., starting to update from p2) The other is a trick employed to use less space by not storing any is_prime data for even numbers)

Let's start with your first two questions. The (2 * i + 3) mapping used in is_prime[i] seems to be a spatial optimization to reduce the space used to half. (i.e., no even number is represented in is_prime list) The mapping helps iterate only the list of odd numbers starting from 3, up to size. If you actually replace the i variable in (2i + 3) with the initial value of size, you will see that you end up with n. (or n-1, depending on whether n is even or odd)

Your third question is relatively more straightforward. In the outer loop, i iterates over the space of odd integers up to n. As there is a mapping of (2i + 3) in is_prime, p is assigned that value. From that point on, p represents the actual prime value which is to be used in the inner loop.

The comment in your fourth question simply further explains the mathematical idea of starting to iterate from p2. As the loop constitutes i to be part of a mapping (to actual values) the p2 is further expressed in terms of that variable i. I think that comment attempts to express the use of 2 * i**2 + 6 * i + 3 to initialize the range of j, but is quite unclear.

To answer your final question, we should consider what j actually represents. j represents the space of odd numbers to be updated. Similar with the loop for i, j iterates not over the actual values, but on the odd numbers. The initial value of j is 2 * i**2 + 6 * i + 3, because when you replace that value with the i variable in (2*i + 3) (i.e., the mapping from the odd numbers' space to the set of actual values), you obtain 4 * i**2 + 12 * i + 9, which is p2.

The inner loop is basically assigning is_prime[j]=False to all the cells that represent the multiples of the actual prime value p, starting from the one representing the value p2.

Answer 2

A naive sieve implementation would have a is_prime array that represents all n numbers we want to check. So its size would be n. Then for each p, we start at 2*p and mark it as "not prime" then go to 3*p, 4*p, 5*p, etc, marking each one as "not prime". So for example, when p = 2, we mark 4, 6, 8, 10, 12, etc as "not prime". Then when p = 3 we mark 6, 9, 12, 15 as "not prime".

I suggest that you implement this algorithm yourself to understand how it works before moving on to implementations. The code you are looking at uses some tricks to reduce the work done. But to understand those tricks, you need to understand the base algorithm.

how did they come up with the formula for the size

This comes from solving the formula n = i * 2 + 3 for i where n is the largest number we will check for primeness. It gives an upper bound for the value of i for all of the numbers we want to test.

how did they get p = i * 2 + 3

This allows us to test only odd numbers starting at 3. Note that even numbers are not prime, so we can easily skip them with this formula.

what does the following mean Sieving from p^2, where p^2 = (4i^2 + 12i + 9). The index in is_prime is (2i^2 + 6i + 3) because is_prime[i] represents 2i + 3

Notice how in our naive algorithm, we marked 6 and 12 as "not prime" twice. We clearly do some extra work here. We can avoid this extra work by realizing that for p, we have already marked all composite numbers less than p^2 as "not prime" when we determined each prime less than p.

So we only have to start at p^2 instead of at p. Now for p = 2, we mark 4, 6, 8, 10, 12, etc as before. But then for p = 3, we mark 9, 12, 15, 18, etc, avoiding the double work of marking 6 as "not prime". For these two examples, the amount of double-marking avoided is pretty small, but as p gets larger, this technique adds up to a significant performance increase.

As for the formula p^2 = (4i^2 + 12i + 9), you can derive this from what we call the FOIL method when multiplying (2*i+3)*(2*i+3). For your code, this doesn't really matter, because if you do p = 2*i + 3, then you can calculate p*p directly without worrying about the underlying algebraic manipulations.

why does the range of j begin with 2 * i**2 + 6 * i + 3

We have p^2 = (4i^2 + 12i + 9) and we need to find the index j in is_prime where p^2 = j * 2 + 3. We set these equal and solve for j.