Understanding Sieve Of Eratosthenes In Python

I've found an example code in python that gives out all prime numbers upto n but I simply don't get it, Why does it does what it does? I've read the wikipedia article about the Sie

Solution 1:

Since no one has yet to show a true sieve or explain it, I will try.

The basic method is to start counting at 2 and eliminate 2*2 and all higher multiples of 2 (ie 4, 6, 8...) since none of them can be prime. 3 survived the first round so it is prime and now we eliminate 3*3 and all higher multiples of 3 (ie 9, 12, 15...). 4 was eliminated, 5 survived etc. The squaring of each prime is an optimization that makes use of the fact that all smaller multiples of each new prime will have been eliminated in previous rounds. Only the prime numbers will be left as you count and eliminate non-primes using this process.

Here is a very simple version, notice it does not use modulo division or roots:

defprimes(n): # Sieve of Eratosthenes
    prime, sieve = [], set()
    for q in xrange(2, n+1):
        if q notin sieve:
            prime.append(q)
            sieve.update(range(q*q, n+1, q))
    return prime

>>> primes(100)
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 7379, 83, 89, 97]

The simple approach above is surprisingly fast but does not make use of the fact that primes can be only odd numbers.

Here is a generator based version that is faster than any other I have found but hits a Python memory limit at n = 10**8 on my machine.

defpgen(n): # Fastest Eratosthenes generatoryield2
    sieve = set()
    for q in xrange(3, n+1, 2):
        if q notin sieve:
            yield q
            sieve.update(range(q*q, n+1, q+q))

>>> timeit('n in pgen(n)', setup="from __main__ import pgen; n=10**6", number=10)
5.987867565927445

Here is a slightly slower but much more memory efficient generator version:

defpgen(maxnum): # Sieve of Eratosthenes generatoryield2
    np_f = {}
    for q in xrange(3, maxnum+1, 2):
        f = np_f.pop(q, None)
        if f:
            while f != np_f.setdefault(q+f, f):
                q += f
        else:
            yield q
            np = q*q
            if np < maxnum:
                np_f[np] = q+q

>>> timeit('n in pgen(n)', setup="from __main__ import pgen; n=10**6", number=10)
7.420101730225724>>> list(pgen(10))
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]

To test if a number is prime just do:

>>> 539in pgen(539)
False>>> 541in pgen(541)
True

Here are some hints as to how this more memory efficient version works. It uses a dict to store only the bare minimum of information, the next non-prime numbers (as keys) along with their factors (as values). As each non-prime is found in the dict, it is removed and the next non-prime key is added with the same factor value.

Solution 2:

Firstly, this is not a sieve.

This is how it works. pp is the number that we are going to test. In each iteration of the while loop, we go over all the known primes (ps) and check if they divide pp. If one of them does, pp is not a prime, and we move to the next number. Otherwise, we add pp to the list of primes before moving on.

The line pp%a==0 is basically saying "the remander of pp when divided by a is zero", ie a divides pp and pp is not prime.

This continues until the number we are checking is larger than some upper limit that we have set (lim)

[EDIT: this is a sieve]

isPrime = [Truefor i inrange(lim)]
isPrime[0] = False
isPrime[1] = Falsefor i inrange(lim):
    if isPrime[i]:
        for n inrange(2*i, lim, i):
            isPrime[n] = False

This is not the most efficient sieve (more efficient ones do things in the for n in range(2*i, lim, i): line, but it will work, and isPrime[i] will be true iff i is prime.

Solution 3:

That code is an attempt at using trial division to produce a sequence of primes.

To correct it:

pp = 2
ps = [pp]
lim = raw_input("Generate prime numbers up to what number? : ")
while pp < int(lim):
    pp += 1for a in ps:
        if pp%a==0:
            breakelse:                # unindent
        ps.append(pp)    #  this

To make it much more efficient (in fact, optimal) trial division:

pp = 2
ps = [pp]
lim = raw_input("Generate prime numbers up to what number? : ")
while pp < int(lim):
    pp += 1for a in ps:
        if a*a > pp:         # stop
            ps.append(pp)    #  earlybreakif pp%a==0:
            break

Solution 4:

The above implementation produces wrong answers. I've done some changes to the code.

But, here's how the above code works.

pp = 2ps = [pp]

We know that the first prime number is 2, so, we generate a list containing only the number 2.

lim = raw_input("Generate prime numbers up to what number? : ")

The above line takes an input from the user, which gives us the upper limit of the prime numbers to generate.

while pp < int(lim):    # 1
      pp += 1# 2
      primeFlag = True# 3for a in ps:      # 4if pp%a==0:
             primeFlag = Falsebreakif primeFlag:     # 5
          ps.append(pp)

The numbered lines do the following things.

1. Runs a loop until the upper limit is reached.
2. Increments the pp variable by 1.
3. Sets a flag variable which is used for testing if the number is prime.
4. This for loop iterates over the list of prime numbers stored in ps and checks that the current number, pp is divisible by any one of those numbers, if yes, then the number is not prime and the primeFlag is set to False and we break out of the inner for loop.
5. If the number was not divisible by any of the primes before it, then it must be a prime, hence, the variable primeFlag is True and the if statement appends the list ps with pp.