Fast computation of prime numbers in Python

Prime Numbers
Before I present the Python code, let's quickly introduce prime numbers. A prime number is a number than can only be divided by itself and 1. For example:

• 7 is a prime number since it can only by divided by either 1 and 7
• 6 is not a prime number since it can divided by 1,2,3,6

The Greeks
Now, let's pay tribute to the Greeks, and for this occasion, one in particular: Eratosthenes (276 BC – c. 195 BC). Eratosthenes is famous for its sieve, which is basically an algorithm to find prime numbers. The idea is very simple for our times but one has to marvel at the genius of having it more than 2000 years ago. His idea was for each number n, mark each number greater than n and dividable by n to be non-prime. So technically what the algorithm does, is for each number p > n, find the remainder of p/n, if the remainder is 0, then p is not a prime. Technically rather than testing every number, Eratosthenes sieve marks directly all the products of a given prime number as non-prime.

Python, what else ?
Here's the code to compute the prime numbers contained in the first 100,000 first natural numbers, it doesn't even take a second to complete. I'm making use of the matrix package Numpy previously introduced in my blog entry about fractals.

import numpy as N

pr = N.ones(100000,'bool') # Interval to search for primes
pr[0] = pr[1] = False
pr[4::2] = False
for x in xrange(3,len(pr),2):
    if pr[x]:           # If x is a prime
        pr[::x] = False # Mark every x entry as 0 (non-prime)
        pr[x]   = True  # Except the one for x

# Finally filter out zero entries
pr2 = pr.nonzero()[0]

And that's all there is to it.

Last but not least, here's a link for verifying big (and small) prime numbers: http://www.numberempire.com/primenumbers.php. This website can verify numbers up to $10^{128}$ which is much superior to any other website I visited.