Sieve Of Eratosthenes In Python 


My latest posts can be found here: Previous blog posts:
Additionally, some earlier writings: 
The Sieve of Eratosthenes

Source: Wikipedia[0] 
So how does the Sieve of Eratosthenes work? We start by writing a flag for every number from 2 up to the limit we're interesting in.
We set these flags as follows:
#!/usr/bin/python limit = 121 flags = {} for i in range(2,limit+1): flags[i] = True p = 2 while p <= limit: if flags[p] == True: print p m = p*p while m <= limit: flags[m] = False m += p p += 1 
Now we move to the next number and go again. If its flag is True then the number is prime, so we print it, and step forward setting the flags for its multiples.
At each stage we know that every composite number up to $p^2$ has definitely been marked as such. So the next unmarked number is always going to be prime, so we call that $p,$ print it, mark all its multiples, and so on.
An important feature of this algorithm is that it doesn't need us to do any division. Until comparatively recently routines called the "Sieve of Eratosthenes" actually used trial division, and so weren't the original algorithm at all. We can see here though that division is not necessary.
The difficulty, however, is that we need to start by initialising a flag for every number up to the limit we're interesting in. If we don't know how far we're going to go, or if we're limited on space, this is a problem. However, we can overcome both these problems by using a more dynamic version.
#!/usr/bin/python from collections import defaultdict # "cmps" is our dictionary # of known composites def create_list(): return [] cmps = defaultdict(create_list) def generate_primes(): v = 2 while True: yield v cmps[v*v].append(v) v += 1 while v in cmps: for f in cmps[v]: cmps[v+f].append(f) del cmps[v] v += 1 
When we look at the next candidate, $v,$ we see if it is in our dictionary of composites. If so then we do not return it, because we know it's not prime. Instead, we look at all the factors $f$ we know of it, and for each we insert into our dictionary that $v+f$ has $f$ as a factor. Having done so we can delete $v$ from our dictionary, thus saving space (and time for subsequent lookups).
We continue until we find a number that is not in our dictionary of composites, and realise that it must be prime. Thus we return it.
def generate_primes(): yield 2 v = 3 while True: yield v cmps[v*v].append(v) v += 2 while v in cmps: for f in cmps[v]: cmps[v+f+f].append(f) del cmps[v] v += 2 
We could go further and also treat "3" as a special case  we can save that for later in case we need it. It's definitely more convoluted to do so, and may not be worth the effort.
n = 2 for p in generate_primes(): while n < p: if is_perrin_q(n): print n n += 1 if clock() >= 100: break n += 1 
When $n$ gets to equal $p$ we don't need to perform a test, because we know that it will pass (because the Perrin test is always true for prime numbers). So we step to the next potential number, and go round the loop with the next prime.
So, timing?

As you can see we've got nearly twice as far, and now we're spending only about 6% of our time in the prime generation routine. So once again it's the Perrin testing that's the slowest part, so once again we ask, how can we make that faster.
That's for next time.
Fast Perrin Test  :  The Point Of The Banach Tarski Theorem ... 
You should follow me on twitter 
I've decided no longer to include comments directly via the Disqus (or any other) system. Instead, I'd be more than delighted to get emails from people who wish to make comments or engage in discussion. Comments will then be integrated into the page as and when they are appropriate.
If the number of emails/comments gets too large to handle then I might return to a semiautomated system. That's looking increasingly unlikely.
Quotation from Tim BernersLee 