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 !