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.

A word about fractals and Python

Introducing Python
There's a nice language that I use pretty much for everything, it is Python. What's so special about Python ? Well it's simple, easy to learn and yet powerful. Its syntax is particularly clear and it can also be compiled into Java bytecode thanks to Jython. I've been using it for years and pretty much anything else I looked at never convinced me as much.

One limitation of Python though is the execution speed, and I'm talking about the pure Python not Jython which is much faster. Since in the next part of this little blog post I am going to do some pretty intensive computations, namely fractals, it's important to mention the speed aspect of things. But, fear not, python is expanded with tons of packages which are C++ fast (which is often as fast as you get, excepting pure machine code which can also be derived from Python but that's another story).

Introducing Numpy and Pylab
Numpy is a mathematic package for python that binds to a huge library of mathematical primitives. Basics include matrices and operations to manipulate them. It also includes more advanced capabilities such as matrix inversion, least square computation, singular value decomposition and a lot more.
Pylab or matplotlib is a superb package to display curves and images. This thing looks good, really good. So check the link if you're looking for some beautiful plotting in a very few lines of code.

Let's enter the fractal world
First let's have a little teaser of a fractal image computed with the program I'm going to present later:
Let's get acquainted with the Mandelbrot set. The Mandelbrot fractal is computed by looking at the divergence of the complex (see this page for an explanation of complex numbers) function $z_{n}=z_{n-1}^{2}+c$ with the initial condition $z_{0}=c$. Mandelbrot's set is computed by computing a z function at every position (x,y) of the 2D space, and the constant c is defined as c=x+j.y.

Next to determine the color of the pixel, we have to determine whether the function z has diverged. We know for sure that it will diverge if $\left|z_{n}\right|\geq4$. So the first n for which the previous condition is true determines the color of the point.

Implementation in Python
As mentioned earlier, for optimal computation speed, we're going to make use of matrices to compute all the points simultaneously.

import numpy as N
import pylab as P

dx = ( .663,.007) # X start position, X width
dy = (-.191,.007) # Y start position, Y width
v  = N.zeros((1024,1024),'uint8')    # The color matrix

# Create the constant matrix and initialize it
c  = N.zeros((1024,1024),'complex')
c[:].real = N.linspace(dy[0],dy[0]+dy[1],
                       c.shape[0])[:,N.newaxis]
c[:].imag = N.linspace(dx[0],dx[0]+dx[1],
                       c.shape[1])[N.newaxis,:]
z  = c.copy()           # The z function is initialized from c
for it in xrange(256):  # Use 256 colors
    z *= z              # Compute z = z*z
    z += c              # Compute z = z + c
    # Set colors for which z has diverged
    v += (N.abs(z) >= 4)*(v == 0)*it
P.imgshow(v)      # Display the image
P.show()

And voila, this will generate a beautiful fractal picture like the one shown above. Note that you can experiment with different start positions and widths for dx and dy. The smaller the width, the more you zoom in, the larger the more you zoom out. Have fun !

Let's start the new year on a positive note: Twitter sucks !

Now that may come as a pretty assertive statement, almost brutal. But let me explain.

The context
Well you have to put yourself in my position for 2 minutes. Here in the silicon valley, but it's probably the same across USA, when I watch the TV, I keep hearing Twitter this, Twitter that. From the cooking shows, to the news, to the UFC fighting organization, you hear about Twitter everywhere, all the time. It's a perpetual brainwashing machine. It's like the internet bubble that exploded in 1999 has re-inflated itself with Web 2.0 !

The initiation
Eventually, be it due to the perpetual brainwashing or to simple curiosity, I finally opened myself a twitter account. First, when I logged in, the twitter website asked me: Tell what you are doing now. Ahem, isn't that a bit too personal a question ? Ok I get it, that's the concept, exhibitionism. Put your life out like you were a popstar, forget intimacy, just let it go. Well, that's where things got difficult for me, I wasn't doing anything special just typing on my computer and trying to find something smart to say about what I was doing which was "trying to find something smart to say about what I was doing", etc...

Seeing the light
So anyway, using Twitter, I really got to realize it's amazing power. It's like an email, but limited to 130 or so characters (I forgot the exact number anyway). Now I understand why everybody talks about it, what a revolution ! Limiting an email to 130 characters, who would have thought of that. It's incredible we had to wait for Web 2.0 to have such revolutionizing concept !

Size doesn't matter
Now it's time to realize, the best quotes from literature are less than 130 characters. Now, with a completely new open mind, I started to see what people were saying, which basically is: "I'm gardening today, it's sunny", "Let's stay in touch", "We announce a new product", etc...

A drop of water in the ocean
All this Twitter experience finally convinced me that people, in the end, just want to feel that someone cares about them, that somehow in the multitude of the human population, someone is genuinely interested in them. However, thanks to Twitter, now I realize that most pretend to be interested in them simply to receive a reciprocal faked interest.

Now join me on Twitter, come read my blog, visit my website.
You won't be disappointed !

And now for something completely different

That's it.
I've done it.
I'm now part of the big internet soup, ready to share my insignificant self !