Sieve of Eratosthenes (Python, arrays)

From LiteratePrograms

The Sieve of Eratosthenes is an algorithm for rapidly locating all the prime numbers in a certain range of integers. It operates by marking as composite all nontrivial multiples of each prime in sequence until only primes remain unmarked. It is most well-suited to implementation using arrays, which are compact and feature random access.

1 simple implementation
2 more complex implementation
3 putting it all together
4 επίλογος

[edit]

simple implementation

We follow the algorithm outlined in the introduction.

<<sieve>>=
  initialize the sieve
  test each possible factor in sequence
          mark nontrivial multiples
  return unmarked primes

First we set up an array of $n 2$ flags, initially all set to 1, indicating primes.

<<initialize the sieve>>=
flags = zeros(n**2) + 1

Next we take each nontrivial factor, from $[2, n)$ ...

<<test each possible factor in sequence>>=
for i in arange(2,n):

... and mark the nontrivial multiples.

we start at $i 2$ because lower multiples will have been tested earlier in the sequence
by using a stride of $i$ all higher multiples are set in the same statement

<<mark nontrivial multiples>>=
flags[i*i::i] = 0

As a result, the remaining unmarked elements are $0$ , $1$ , and the primes; we drop the first two to return only the primes.

<<return unmarked primes>>=
return nonzero(flags)[2:]

[edit]

more complex implementation

We can attempt to speed the algorithm up by precalculating the first few passes. With a variant of wheel factorization, rather than initializing the sieve to all ones, we replicate a pattern without multiples of 2 and 3, and then reset the flags for the trivial multiples themselves.

<<initialize precalculated sieve>>=
flags = resize((0,1,0,0,0,1), (n**2,))
put(flags, (0,2,3), 1)

Because of the precalculation, we start testing from the next remaining prime, $5$ . Also, instead of testing each factor, prime or composite, we check the results calculated so far to test only the prime factors.

<<test each prime factor in sequence>>=
for i in arange(5,n,2):
    if flags[i]:

The remainder of the algorithm remains unchanged.

<<koskinon>>=
  initialize precalculated sieve
  test each prime factor in sequence
          mark nontrivial multiples
  return unmarked primes

[edit]

putting it all together

Finally, we generate a list of the primes less than 9'999, and check that sieve and koskinon yield identical results.

<<sieve.py>>=
from Numeric import *

def sieve(n):
sieve

def koskinon(n):
koskinon

if __name__ == "__main__":
        print sieve(100)
        print multiply.reduce(koskinon(100)==sieve(100))

which we expect to produce:

[   2    3    5    7   11   13   17   19   23   29   31   37   41   43   47
      53   59   61   67   71   73   79   83   89   97  101  103  107  109
      ...
      9781 9787 9791 9803 9811 9817 9829 9833 9839 9851 9857 9859 9871 9883
      9887 9901 9907 9923 9929 9931 9941 9949 9967 9973]
1

[edit]