Take the 2-minute tour ×
Code Review Stack Exchange is a question and answer site for peer programmer code reviews. It's 100% free, no registration required.

So I've been doing a lot of Project Euler lately and just wanted to make sure my implementation was as good as it could be. Does anyone have any suggestions to speed this up?

def sieve(upperlimit):
    # mark off all multiples of 2 so we can use 2*p as the step for the inner loop
    l = [2] + [x if x % 2 != 0 else 0 for x in range(3, upperlimit + 1)]

    for p in l:
        if p ** 2 > upperlimit:
            break
        elif p:
            for i in range(p * p, upperlimit + 1, 2 * p):
                l[i - 2] = 0
    # filter out non primes from the list, not really that important i could work with a list full of zeros as well
    return [x for x in l if x]
share|improve this question
    
For most project euler problems this version of the sieve together with an implementation of the Miller-Rabin primality test should be all you need. I'm not sure saving 2 places in the array is worth the number of times your going to adjust the index with l[i - 2] though. –  Jackson Feb 21 at 16:11
add comment

4 Answers

Your first setting of all the even numbers to 0 is not very efficient, the whole point of sieving is to avoid those costly modulo operations. Try the following:

l = range(2, upperlimit+1) # use list(range(...)) in Python 3
l[2::2] = [0] * ((len(l) - 3) // 2 + 1)

You can do a similar thing for the setting of zeros of the sieve for other prime numbers, but it gets complicated to figure out how many zeros to add on the right.

share|improve this answer
    
Oh nice catch, yea that saved me a decent chunk of time when generating the list –  Igglyboo Feb 21 at 16:27
1  
Note that this fails in python3 (and the question is tagged python3). –  Bakuriu Feb 21 at 17:22
    
Yep, Python 3 requires wrapping the range in a list for it to work, added it as a comment above. –  Jaime Feb 21 at 17:59
add comment
  1. Here's my starting point for computing the performance improvement due to the various revisions below: how long does it take to sieve for the prime numbers below 108?

    >>> from timeit import timeit
    >>> test = lambda f: timeit(lambda:f(10**8), number=1)
    >>> t1 = test(sieve)
    

    The exact number is going to depend on how fast your computer is, so I'm going to compute performance ratios, but for the record, here it is:

    >>> t1
    78.9875438772142
    
  2. Your initialization of the list l takes more than half the time, so let's try a cheaper approach. Let's also give this array a better name, and make it a Boolean array while we're about it.

    def sieve2(n):
        """Return a list of the primes below n."""
        prime = [True] * n
        for p in range(3, n, 2):
            if p ** 2 > n:
                break
            if prime[p]:
                for i in range(p * p, n, 2 * p):
                    prime[i] = False
        return [2] + [p for p in range(3, n, 2) if prime[p]]
    

    When optimizing a function like this, it's always worth keeping the un-optimized version around to check the correctness of the optimized version:

    >>> sieve(10**6) == sieve2(10**6)
    True
    

    This already runs in less than a third of the time:

    >>> test(sieve2) / t1
    0.30390444573149544
    
  3. We could avoid the test for p ** 2 > n by computing a tighter limit for the loop. Note that I've used n ** .5 here as this is slightly faster than math.sqrt(n).

    def sieve3(n):
        """Return a list of the primes below n."""
        prime = [False, False, True] + [True, False] * (n // 2)
        for p in range(3, int(n ** .5) + 1, 2):
            if prime[p]:
                for i in range(p * p, n, 2 * p):
                    prime[i] = False
        return [p for p in range(2, n) if prime[p]]
    

    This makes little difference to the overall runtime:

    >>> test(sieve3) / t1
    0.2971086436068156
    
  4. We can accumulate the result as we go, instead of in a separate iteration at the end. Note that I've cached result.append in a local variable to avoid looking it up each time round the loop.

    def sieve4(n):
        """Return a list of the primes below n."""
        prime = [False, False, True] + [True, False] * (n // 2)
        result = [2]
        append = result.append
        sqrt_n = (int(n ** .5) + 1) | 1    # ensure it's odd
        for p in range(3, sqrt_n, 2):
            if prime[p]:
                append(p)
                for i in range(p * p, n, 2 * p):
                    prime[i] = False
        for p in range(sqrt_n, n, 2):
            if prime[p]:
                append(p)
        return result
    

    Again, this makes very little difference:

    >>> test(sieve4) / t1
    0.286016401170129
    
  5. We can use Python's slice assignment instead of a loop when setting the sieve entries to False. This looks wasteful since we create a large list and then throw it away, but this avoids an expensive for loop and the associated Python interpreter overhead.

    def sieve5(n):
        """Return a list of the primes below n."""
        prime = [True] * n
        result = [2]
        append = result.append
        sqrt_n = (int(n ** .5) + 1) | 1    # ensure it's odd
        for p in range(3, sqrt_n, 2):
            if prime[p]:
                append(p)
                prime[p*p::2*p] = [False] * ((n - p*p - 1) // (2*p) + 1)
        for p in range(sqrt_n, n, 2):
            if prime[p]:
                append(p)
        return result
    

    This gives a small but noticeable improvement:

    >>> test(sieve5) / t1
    0.2617646381557855
    
  6. For big improvements to the performance of numerical code, we can use NumPy.

    import numpy
    
    def sieve6(n):
        """Return an array of the primes below n."""
        prime = numpy.ones(n, dtype=numpy.bool)
        prime[:2] = False
        prime[4::2] = False
        sqrt_n = int(n ** .5) + 1
        for p in range(3, sqrt_n, 2):
            if prime[p]:
                prime[p*p::2*p] = False
        return prime.nonzero()[0]
    

    This is more than 6 times as fast as sieve5, and more than 25 times as fast as your original code:

    >>> test(sieve6) / t1
    0.03726392181902129
    
  7. We could avoid allocating space for the even numbers, improving memory locality:

    def sieve7(n):
        """Return an array of the primes below n."""
        prime = numpy.ones(n // 2, dtype=numpy.bool)
        sqrt_n = int(n ** .5) + 1
        for p in range(3, sqrt_n, 2):
            if prime[p // 2]:
                prime[p*p // 2::p] = False
        result = 2 * prime.nonzero()[0] + 1
        result[0] = 2
        return result
    
    >>> test(sieve7) / t1
    0.029220096670965198
    
  8. And finally, an implementation that sieves separately for primes of the form 6𝑖 − 1 and 6𝑖 + 1, due to Robert William Hanks:

    def sieve8(n):
        """Return an array of the primes below n."""
        prime = numpy.ones(n//3 + (n%6==2), dtype=numpy.bool)
        for i in range(3, int(n**.5) + 1, 3):
            if prime[i // 3]:
                p = (i + 1) | 1
                prime[       p*p//3     ::2*p] = False
                prime[p*(p-2*(i&1)+4)//3::2*p] = False
        result = (3 * prime.nonzero()[0] + 1) | 1
        result[0] = 3
        return numpy.r_[2,result]
    

    This is about 40 times as fast as the original sieve:

    >>> test(sieve8) / t1
    0.023447068662434022
    
share|improve this answer
    
I'm curious, if n ** 0.5 is faster than math.sqrt(n), why isn't sqrt() implemented using n ** 0.5? Does math.sqrt(n) have better accuracy or something? –  200_success Feb 21 at 20:59
    
@200_success: sqrt(n) involves a function call and n ** .5 doesn't. –  Gareth Rees Feb 21 at 21:00
    
One function call, not in a loop, isn't really worth fussing about. –  200_success Feb 21 at 21:05
    
7 can be generalized using wheel factorization, which gets you another large improvement. Filling the sieve in chunks that fit into L1 or L2 cache gets another factor 5 or more according to my experience –  Niklas B. Feb 22 at 0:14
    
@NiklasB.: Sounds interesting! Can you post your code? –  Gareth Rees Feb 22 at 17:09
show 1 more comment

You have correctly implemented the optimization p ** 2 > upperlimit, but, there is a more efficient way to do it...

p ** 2 > upperlimit calculates the square of p for every iteration through the loop. The product is not hugely expensive, but it is also totally unnecessary....

You can calculate the square-root of the upperlimit just once, and then reuse that calculated value as a direct comparison for p. Consider:

rootlimit = math.sqrt(upperlimit)
for p in l:
    if p > rootlimit:
        break;
    ....

Additionally, it is a small thing, but if you have break, continue, or return, inside a conditional inside a loop, then there is no need to use elif... It is just a small thing, but your code could be:

rootlimit = math.sqrt(upperlimit)
for p in l:
    if p > rootlimit:
        break;
    if p:
        ....
share|improve this answer
add comment
for p in l:
    if p ** 2 > upperlimit:
        break
    elif p:
        ...

In this loop p is often zero, while p ** 2 > upperlimit is always false until the break. That means you evaluate both conditions always. If you put if p: first, only one condition is evaluated whenever p is not prime.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.