How to speed up Sieve of Eratosthenes python list generator

Question

My problem comes directly from the CS circles site. It's the last problem on the bottom of this page called 'Primed for Takeoff'. The basic rundown is they want a list of 1,000,001 length, where the index of each item is True if the index is prime, while the index of each item is False if it is not prime.

So, for example, isPrime[13] is True. isPrime[14] is False.

The first little bit of the list 'isPrime' would look like:

isPrime = [False, False, True, True, False, True, False, True, False, False, ...]

My problem is the 7-second time limit they have. I have a working code below with a lower number, 101, for debugging purposes. When I bump it up to their required 1,000,001 list length, it takes way too long, I ended up killing the program after a few minutes. This is my working (at length 101), very slow code:

n = 101
c = 2
isPrime = [False,False]
for i in range(2,n):
    isPrime.append(i)

def ifInt(isPrime):
    for item in isPrime:
        if type(item) == int:
            return item
for d in range(2,n):
    if c != None:
        for i in range(c,n,c):
            isPrime[i] = False
        isPrime[c] = True
        c = ifInt(isPrime)

Then there's this one I found here. It has a quicker run time, but only outputs a list of primes, not a list of n length with list[n] returning True for primes and false otherwise.

I'm not sure if this 2nd bit of code really holds the key to creating a prime list of length 1,000,001 in under 7 seconds, but it was the most relevant thing I could find in my research.

def primes_sieve1(limit):
limitn = limit+1
primes = dict()
for i in range(2, limitn): primes[i] = True

for i in primes:
    factors = range(i,limitn, i)
    for f in factors[1:]:
        primes[f] = False
return [i for i in primes if primes[i]==True]

print primes_sieve1(101)

CS circles seems pretty commonly used, but I wasn't able to find a working version of this for Python. Hopefully this will be an easy solve for someone.

This question differs from this one because am not trying to only create a list of primes quickly, but to create a list of all positive integers from 0 to n that are marked as prime by True and non prime by False.

Answer 1

I realized that there a lot of optimizations on SO, but they are rarely ever explained by others for the prime sieve algorithm, so it makes them difficult to approach by beginners or first time creators of the algorithm. All the solutions here are in python, to be on the same page for speed and optimizations. These solutions will progressively become faster and more complex. :)

Vanilla Solution

def primesVanilla(n):
    r = [True] * n
    r[0] = r[1] = False
    for i in xrange(n):
        if r[i]:
            for j in xrange(i+i, n, i):
                r[j] = False
    return r

This is a very straightforward implementation of the Sieve. Please make sure you understand what is going on above before proceeding. The only slight thing to note is that you start marking not-primes at i+i instead of i, but this is rather obvious. (Since you assume that i itself is a prime). To make tests fair, all numbers will be for the list up to 25 million.

real    0m7.663s  
user    0m7.624s  
sys     0m0.036s  

Minor Improvement 1 (Square roots):

I'll try to sort them in terms of straight-forward to less straight-forward changes. Observe that we do not need to iterate to n, but rather only need to go up to the square root of n. The reason being that any composite number under n, must have a prime factor under or equal to the square root of n. When you sieve by hand, you'll notice that all the "unsieved" numbers over the square root of n are by default primes.

Another remark is that you have to be a little bit careful for when the square root turns out to be an integer, so you should add one in this case so it covers it. IE, at n=49, you want to loop until 7 inclusive, or you might conclude that 49 is prime.

def primes1(n):
    r = [True] * n
    r[0] = r[1] = False
    for i in xrange(int(n**0.5+1)):
        if r[i]:
            for j in xrange(i+i, n, i):
                r[j] = False
    return r

real    0m4.615s
user    0m4.572s
sys     0m0.040s

Note that it's quite a bit faster. When you think about it, you're looping only until the square root, so what would take 25 million top level iterations now is only 5000 top level.

Minor Improvement 2 (Skipping in inner loop):

Observe that in the inner loop, instead of starting from i+i, we can start from i*i. This follows from a very similar argument to the square root thing, but the big idea is that any composites between i and i*i have already been marked by smaller primes.

def primes2(n):
    r = [True] * n
    r[0] = r[1] = False
    for i in xrange(int(n**0.5+1)):
        if r[i]:
            for j in xrange(i*i, n, i):
                r[j]=False
    return r

real    0m4.559s
user    0m4.500s
sys     0m0.056s

Well that's a bit disappointing. But hey, it's still faster.

Somewhat Major Improvement 3 (Even skips):

The idea here is that we can premark all the even indices, and then skip iterations by 2 in the main loop. After that we can start the outer loop at 3, and the inner loop can skip by 2*i instead. (Since going by i instead implies that it'll be even, (i+i) (i+i+i+i) etc.)

def primes3(n):
    r = [True] * n
    r[0] = r[1] = False
    for i in xrange(4,n,2):
        r[i] = False    
    for i in xrange(3, int(n**0.5+1), 2):
        if r[i]:
            for j in xrange(i*i, n, 2*i):
                r[j] = False
    return r

real    0m2.916s
user    0m2.872s
sys     0m0.040s

Cool Improvements 4 (wim's idea):

This solution is a rather advanced trick. Slice assignment is faster than looping, so this uses python's slice notation: r[begin:end:skip]

def primes4(n):
    r = [True] * n
    r[0] = r[1] = False 
    r[4::2] = [False] * len(r[4::2])
    for i in xrange(3, int(1 + n**0.5), 2):
        if r[i]:
            r[i*i::2*i] = [False] * len(r[i*i::2*i])
    return r

10 loops, best of 3: 1.1 sec per loop

Slight Improvement 5

Note that python reslices the r[4::2] when it calculates the length, so this takes quite a bit of extra time since all we need for it is computing the length. We do use some nasty math to achieve this though.

def primes5(n):
    r = [True] * n
    r[0] = r[1] = False 
    r[4::2] = [False] * ((n+1)/2-2)
    for i in xrange(3, int(1 + n**0.5), 2):
        if r[i]:
            r[i*i::2*i] = [False] * ((n+2*i-1-i*i)/(2*i))
    return r

10 loops, best of 3: 767 msec per loop

Assignment speed-up (Padraic Cunningham):

Note that we assign an array with all True and then set half (the evens) to be False. We can actually just start with a boolean array that's alternating.

def primes6(n):
    r = [False, True] * (n//2) + [True]
    r[1], r[2] = False, True
    for i in xrange(3, int(1 + n**0.5), 2):
        if r[i]:
            r[i*i::2*i] = [False] * ((n+2*i-1-i*i)/(2*i))
    return r

10 loops, best of 3: 717 msec per loop

Don't quote me on this, but I think without some nasty math methods, there is no obvious improvements to this last version. One cute property that I tried, but did not turn out to be any faster is noting that primes other than 2,3 must be of the form 6k+1 or 6k-1. (Note that if it's 6k, then divisible by 6, 6k+2 | 2, 6k+3 | 3, 6k+ 4 | 2, 6k+5 is congruent to -1 mod 6. This suggests that we can skip by 6 each time and check both sides. Either from a poor implementation on my side, or python internals, I was unable to find any meaningful speed increase. :(