Prime Numbers in Python

Here is the full generator. The source can also be downloaded from here: primegen.py
It's based on the code posted here

Output:

Performance

Our performance data has been generated from profiling.zip using pypy and python2. primegen.py will need to be included in the same folder. The generators used for profiling were taken from the page posted above, and each generator is denoted by the last name of its author.

The following table shows the time and memory used to iterate over all primes less than 10^8. The ability to skip more composites per iteration, and inlining adding factors to comp, allows our algorithm to perform about 36% faster than the next best.

The following graph shows the time taken to iterate up to 10^7. Even though these algorithms involve finding a factor of n, their performance appears linear with respect to n.

Explanation

The algorithm relies on having a rolling dictionary of upcoming composite numbers, and being able to retrieve their factors quickly. We'll first analyze a simple version of this dictionary and then look at ways to speed it up.

Our dictionary of composites, hence named comp, will return a factor of n if n is composite, otherwise it will have no entry for n. That is, if n is composite, comp[n]=f for n%f=0, otherwise n is prime and comp[n] will be missing.

We will build comp as we iterate over n=2,3,4,... For every iteration we will find a factor, f, of n using f=comp[n]. If n is prime, we'll take f=n. Otherwise, the factor will be used to calculate the next composite not already in comp that's a multiple of f. That is, we'll iterate over values of n+k*f for k=1,2,3,... until we find a value not in comp. This will give us the next entry comp[n+k*d]=f.

The code for this is much more succinct than the derivation:

Indeed, this version is sufficient for a simple prime number generator:

The first big improvement comes from Will Ness, who noted that prime factors are being added to comp before they need to be. That is, for a prime p, we don't need to be tracking p as a factor of n until n=p^2. For n<p^2, there will be some other prime factor of n that will identify n as a composite.

For example, consider p=7. We can see that for n=8,9,...,48, all values will either be prime, or have 2, 3, or 5 as a factor:

It isn't until n=49 that the smallest factor is 7. Using the simple generator we have so far, let us see what comp is when n=49:

There are a total of 15 entries. We know that we only need to have 4 at most: 2, 3, 5, and 7. These extra entries increase the values skipped over in the n+k*f loop and increase the number of entries that need to be ignored when we call if n in comp.

We can improve on this by delaying the addition of p until we pass n>=p^2, and just rely on the other factors to identify composites. The main problem of this is tracking what the next p will be when n=p^2, which is the exact problem we're already trying to solve with our prime generator!

Fortunately, we can solve this easily by making a recursive call to the primegen iterator we already have. This new iterator will solely be used for tracking the next factor we need to add. We need to be careful about making an infinite number of recursive calls to primegen however, so we're going to manually output the first 2 primes first and then use our main loop.

Here, sq=r^2 where sq is the square we're watching for and r is the next prime we'll add to comp. If we check what comp contains when n=49, it's what we'd expect: {50: 5, 51: 3, 52: 2}. Note that 7 is not an entry because it is implicitly stored in r.

The last remaining way to speed up the algorithm is skip numbers that we know have a factor. For instance, every even number greater than 2 is a composite divisible by 2. Every third number is divisible by 3, etc... In general, if we choose to skip composites divisible by primes 2,3,...,p, then we will be able to skip a larger and larger fraction of integers described by the formula: y=1-(1/2)*(2/3)*(4/5)*...*((p-1)/p).

For the purpose of this generator, we choose to skip integers divisible by 2, 3, and 5. By the table above, this will allow us to skip 73% of all integers we need to test. To do this we'll create a table that tells us how far to jump from n so that n+jump is not divisible by 2, 3, or 5. To simplify this table, note that r is divisible by 2, 3, or 5, if and only if r+30 is also divisible by 2, 3, or 5. Thus, to calculate how far from n to jump so n+jump isn't divisible by 2, 3, or 5, we may write n as n=q*30+r for some 0<r<=30 and use r to find the jump value.

To derive the jump table, we begin by looking at the values 0 to 29:

We now cross off all values divisible by 2, 3, or 5, leaving:

Our jump table will be calculated as follows: given some entry on the table, how many spaces do we have to jump to reach the next value not divisible by 2, 3, or 5. We have

Or, in python:

We can now calculate the next value of n with n+=jump[n%30].

Before we simply drop this table into the main loop of our prime generator, we need to modify the section that adds upcoming composites to comp. Previously, when we added values to comp, we knew that n would iterate over all positive integers without skipping any values. Because of our new jump table, we introduce the possibility that we'll jump over a composite that comp is expecting us to hit.

To fix this section, we'll make use of the jump table. Our fixed loop must jump over all the values that n will jump over, and still be a multiple of f. We can do this by factoring f out of n, and then jumping over values when we check for composites already in comp. Our new loop will look like so:

We'll present what our generator currently looks like, but note that it does not work correctly.

The last modification we need to make is to the primes we manually yield. We must yield at least 2, 3, and 5 to prevent them from being skipped by jump, and we must also yield 7 to prevent primegen() from making infinite recursive calls. We also want to yield as many as possible to minimize the number of recursive calls, but we can't yield any prime greater than 7^2=49. Thus, we'll yield all primes up 49. We'll also rearrange some of the code and add convenience checks for the maximum value of n and maximum number of primes to generate.

Here is the final version of our prime generator:

Notes

We have explained how the prime generating algorithm works and shown that it is faster than current related algorithms. We now list some current ideas for improving it.

One trivial improvement is to increase the primes we jump around (2,3, and 5). This has the drawback of increasing the jump table size for every prime we add, and every additional prime skips fewer and fewer numbers.

When we are querying if the current value of n is composite, f=comp.pop(n,0), n will either be the minimum key of comp or it will not be in comp at all. Any data structure which could exploit this will allow for a nice speedup.

We could potentially speed up the algorithm by precalculating the jump table used when calculating the next composite of f. That is

In practice, the table precalculation adds more time than is saves (6.88s vs the current 4.65s).

Another algorithm we could use relies on primality testing. It would look like so:

It is conceptually simple, uses constant memory, and takes O(log(n)^2) time using the Rabin-Miller primality test. In practice, the isprime() function is too slow to compete with the current primegen for real-world values of n.