Python Eratosthenes Sieve Algorithm Optimization

Question

I'm attempting to implement the Sieve of Eratosthenes. The output seems to be correct (minus "2" that needs to be added) but if the input to the function is larger than 100k or so it seems to take an inordinate amount of time. What are ways that I can optimize this function?

def sieveErato(n):
     numberList = range(3,n,2)

     for item in range(int(math.sqrt(len(numberList)))):
            divisor = numberList[item]
            for thing in numberList:
                    if(thing % divisor == 0) and thing != divisor:
                            numberList.remove(thing)
    return numberList

Answer 1

Your algorithm is not the Sieve of Eratosthenes. You perform trial division (the modulus operator) instead of crossing-off multiples, as Eratosthenes did over two thousand years ago. Here is an explanation of the true sieving algorithm, and shown below is my simple, straight forward implementation, which returns a list of primes not exceeding n:

def sieve(n):
    m = (n-1) // 2
    b = [True]*m
    i,p,ps = 0,3,[2]
    while p*p < n:
        if b[i]:
            ps.append(p)
            j = 2*i*i + 6*i + 3
            while j < m:
                b[j] = False
                j = j + 2*i + 3
        i+=1; p+=2
    while i < m:
        if b[i]:
            ps.append(p)
        i+=1; p+=2
    return ps

We sieve only on the odd numbers, stopping at the square root of n. The odd-looking calculations on j map between the integers being sieved 3, 5, 7, 9, ... and indexes 0, 1, 2, 3, ... in the b array of bits.

You can see this function in action at http://ideone.com/YTaMB, where it computes the primes to a million in less than a second.

Answer 2

You can try the same way Eratosthenes did. Take an array with all numbers you need to check order ascending, go to number 2 and mark it. Now scratch every second number till the end of the array. Then go to 3 and mark it. After that scratch every third number . Then go to 4 - it is already scratched, so skip it. Repeat this for every n+1 which is not already scratched.

In the end, the marked numbers are the prime one. This algorithm is faster, but sometimes need lots of memory. You can optimize it a little by drop all even numbers (cause they are not prime) and add 2 manually to the list. This will twist the logic a little, but will take half the memory.

Here is an illustration of what I'm talking about: http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes

Answer 3

Warning: removing elements from an iterator while iterating on it can be dengerous...

You could make the

    if(thing % divisor == 0) and thing != divisor:

test lighter by splitting it in the loop that breaks when you arrive to the index of 'divisor' and then the test:

for thing in numberList_fromDivisorOn:
    if(thing % divisor == 0):
        numberList.remove(thing)

Answer 4

I followed this link: Sieve of Eratosthenes - Finding Primes Python as suggested by @MAK and I've found that the accepted answer could be improved with an idea I've found in your code:

def primes_sieve2(limit):
    a = [True] * limit               # Initialize the primality list
    a[0] = a[1] = False
    sqrt = int(math.sqrt(limit))+1
    for i in xrange(sqrt):
        isprime = a[i]
        if isprime:
            yield i
            for n in xrange(i*i, limit, i):     # Mark factors non-prime
                a[n] = False
    for (i, isprime) in enumerate(a[sqrt:]):
        if isprime:
            yield i+sqrt

Answer 5

if given unlimited memory and time, the following code will print all the prime numbers. and it'll do it without using trial division. it is based on the haskell code in the paper: The Genuine Sieve of Eratosthenes by Melissa E. O'Neill

from heapq import heappush, heappop, heapreplace
def sieve():
    w = [2,4,2,4,6,2,6,4,2,4,6,6,2,6,4,2,6,4,6,8,4,2,4,2,4,8,6,4,6,2,4,6,2,6,6,4,2,4,6,2,6,4,2,4,2,10,2,10]
    for p in [2,3,5,7]: print p
    n,o = 11,0
    t = []
    l = len(w)
    p = n
    heappush(t, (p*p, n,o,p))
    print p
    while True:
        n,o = n+w[o],(o+1)%l
        p = n
        if not t[0][0] <= p:
            heappush(t, (p*p, n,o,p))
            print p
            continue
        while t[0][0] <= p:
            _, b,c,d = t[0]
            b,c = b+w[c],(c+1)%l
            heapreplace(t, (b*d, b,c,d))
sieve()

Answer 6

This code takes 2 seconds to generate primes less than 10M (it is not mine, i found it somewer on google)

def erat_sieve(bound):
    if bound < 2:
        return []
    max_ndx = (bound - 1) // 2
    sieve = [True] * (max_ndx + 1)
    #loop up to square root
    for ndx in range(int(bound ** 0.5) // 2):
        # check for prime
        if sieve[ndx]:
            # unmark all odd multiples of the prime
            num = ndx * 2 + 3
            sieve[ndx+num:max_ndx:num] = [False] * ((max_ndx-ndx-num-1)//num + 1)
    # translate into numbers
    return [2] + [ndx * 2 + 3 for ndx in range(max_ndx) if sieve[ndx]]