I found a good homework and posted at Codenone since last week. The first question is about Mandelbrot Set. Actually, this problem could be categorized in compute intensive which should be solved in low-level language to let it run as fast as possible. However, we are talking about alternative languages aka scripting languages like Python and Ruby so I will try to solve this problem in Python. For the solution in Ruby, please go to another post.

Since performance issue is discussed in Ruby version, I will try my best to do the same thing in Python. First of all, the main loop of the benchmark is as follow.

def bm(is_mandel):
    step = 0.005
    x = 0.0
    while x < 1.0:
        y = 0.0
        while y < 1.0:
            is_mandel(x,y)
            y += step
        x += step

The bm() will call is_mandel() 40,000 times (200x200). It was converted directly from the Ruby version. Then I implemented the straightforward version using built-in complex type.

def is_mandel_1(x,y,maxiter=30,limit=2):
    p = complex(x,y)
    i = 0
    z = 0+0j
    while abs(z) < limit and i < maxiter:
        z = z*z+p
        i += 1
    return i

Then I used the optimized code from wikipedia so below code don't use complex anymore.

def is_mandel_2(x,y,maxiter=30,limit=2):
    x0,y0 = x,y
    x2,y2 = x*x,y*y
 
    i = 0
    while x2+y2 < limit*limit and i < maxiter:
        x,y = x2-y2+x0,2*x*y+y0
        x2,y2 = x*x,y*y
        i += 1
    return i

There are some overhead regarding limit*limit since Python seems to never optimize this constant. The third is just to calculate lim = limit*limit.

def is_mandel_3(x,y,maxiter=30,limit=2):
    x0,y0 = x,y
    x2,y2 = x*x,y*y
    lim = limit*limit
 
    i = 0
    while x2+y2 < lim and i < maxiter:
        x,y = x2-y2+x0,2*x*y+y0
        x2,y2 = x*x,y*y
        i += 1
    return i

The last try is to implement is_mandel() in C. Sadly, the direct port of inline module like PyInline seems to be freezed since 3 years ago and I don't find PyInline in Ubuntu Feisty. However, there is another alternative module called Pyrex. The syntax is very similar to Python but includes variable type in the code. As a result, I tried to implement is_mandel_4() in Pyrex as follow, says mandelbrot.pyx.

def is_mandel_4(double x,double y,int maxiter=30,int limit=2):
    cdef double x0,y0,x2,y2,lim
 
    x0,y0 = x,y
    x2,y2 = x*x,y*y
    lim = limit*limit
 
    i = 0
    while x2+y2 < lim and i < maxiter:
        x,y = x2-y2+x0,2*x*y+y0
        x2,y2 = x*x,y*y
        i = i+1
    return i

Pyrex will generate mandelbrot.c for later compiling by native compiler like gcc so I wrote a Makefile as below.

PYINCLUDE = -I/usr/include/python$(shell python -c "import sys;print sys.version.split()[0]")
 
all:    mandelbrot.so
 
%.c:    %.pyx
        pyrexc $<
 
%.o:    %.c
        gcc -c -fPIC $(PYINCLUDE) $<
 
%.so:   %.o
        gcc -shared $< -lm -o $@
 
clean:
        @rm -f *.c *.o *.so *~ core core.*

Before using the code, I had to generate and compile by make all. After that I got mandelbrot.so and the code to run that function looks like below.

import mandelbrot
is_mandel_4 = mandelbrot.is_mandel_4

Now it is time to see the truth. I ran above functions in IBM ThinkPad T43 which is Pentium M 1.86 MHz and got below result.

              real   user    sys
is_mandel_1 0.671s 0.616s 0.036s
is_mandel_2 0.879s 0.852s 0.008s
is_mandel_3 0.852s 0.808s 0.020s
is_mandel_4 0.225s 0.180s 0.032s

Surprisingly, complex is faster than optimized float. Other results are similar to my expectation.

Tags: codenone, python, ruby, mandelbrot set, pyinline, pyrex