OXFORD UNIVERSITY  COMPUTING LABORATORY

Fraskell

Pretty pictures

To the untrained eye Fraskell seems to be just a Haskell program that generates images of the Mandelbrot set. Compile it with ye olde Haskell 98 compiler (give or take a flag to enable cpp and, for the graphical front-end, Gtk+HS (for which you may also need this patch)) and you will be able to explore the world of the Mandelbrot set, as shown below, or produces images such as this larger image (1152x864 pixels) from the command line.

Screenshot of Fraskell's GTK Mandelbrot exploring interface

Nothing hugely new or exciting there...that's coming later. But first a brief interlude about the Mandelbrot set to prepare us for what is to come.

The Mandelbrot set

The Mandelbrot set was discovered in 1980 by Benoit B. Mandelbrot. The underlying mathematics are extremely simple - for each point (x, y) on a plane:

  • Let z = 0
  • Repeatedly let z = z*z + (x + y*i)

If the absolute value of z ever becomes greater than 2 then the point (x, y) is not in the Mandelbrot set. Otherwise it is. Simple as that!

The bound of 2 has not been chosen at random. A few moments thought by those with a mathematical background will be sufficient to convince them that once the absolute value of z reaches a value greater than 2 it will stay above 2 forever. A corollary of this is that all points in the Mandelbrot set are within the circle of radius 2 centered on the origin (inclusive).

If set membership were all that was done then monochrome pictures would be all that we could produce, yet pictures of the Mandelbrot set tend to be very colourful. The colours represent, for those points not in the set, how many iterations of the second step above were needed to get a value of z with magnitude greater than 2. In the picture above the more blue a point the greater the number of iterations were needed.

Of course it is not feasible for real implementations to continue iterating until set membership has been determined - it could take an enormous number of iterations requiring far too much CPU time to be feasible, and indeed has been conjectured to be undecidable. Implementations thus tend to try iterating at most only a fixed number of times, assuming the point to be in the set if it hasn't been shown not to be by that point.

Speed

So now we know how the Mandelbrot set is defined; we have a program that draws pretty pictures of it for us; only one thing is left to do - make it faster!

The obvious way, in a functional language, to implement the inner core of a Mandelbrot program is with a simple recursive function that passes the iteration count so far to recursive calls and returns a value indicating what colour a point should be coloured. For example, consider this simplified version, where the colour is 0, 1 or 2 for 1, 2 and 3 iterations respectively and 0 if more iterations than that are needed:

mb n z (x, y) = if n > 2
                then 0
                else let z' = z*z + x + y*i
                     in  if |z'| > 2
                         then n
                         else mb (n + 1) z' (x, y) 

Now note that we always call this function with n equal to 0 (except in the recursive calls). We could therefore equally well call mb0 below, which is the same as mb except specialised for the case when n equals 0. Note that the outer if statement has completely gone as we know that the first branch can never be chosen.

mb0 z (x, y) = let z' = z*z + x + y*i
               in  if |z'| > 2
                    then 0
                    else mb 1 z' (x, y) 

But we can perform the equivalent specialisation for the case where n is 1 and substitute the code where the call to mb is made (note the alpha renaming and variable capture):

mb0 z (x, y) = let z' = z*z + x + y*i
               in  if |z'| > 2
                    then 0
                    else if 1 > 2
                         then 0
                         else let z'' = z'*z' + x + y*i
                              in  if |z'| > 2
                                  then 1
                                  else mb (1 + 1) z' (x, y)

Again we know the else clause of the outer if statement of the new code will always be chosen so we can simplify it away:

mb0 z (x, y) = let z' = z*z + x + y*i
               in  if |z'| > 2
                    then 0
                    else let z'' = z'*z' + x + y*i
                         in  if |z'| > 2
                             then 1
                             else mb 2 z' (x, y)

Again we can replace the call to mb with a specialised textual copy and simplify:

mb0 z (x, y) = let z' = z*z + x + y*i
               in  if |z'| > 2
                    then 0
                    else let z'' = z'*z' + x + y*i
                         in  if |z''| > 2
                             then 1
                             else let z''' = z''*z'' + x + y*i
                                  in  if |z'''| > 2
                                  then 2
                                  else mb 3 z' (x, y)

This time after we substitute the code for mb we get this:

mb0 z (x, y) = let z' = z*z + x + y*i
               in  if |z'| > 2
                    then 0
                    else let z'' = z'*z' + x + y*i
                         in  if |z''| > 2
                             then 1
                             else let z''' = z''*z'' + x + y*i
                                  in  if |z'''| > 2
                                  then 2
                                  else if 3 > 2
                                       then 0
                                       else let z' = z*z + x + y*i
                                            in  if |z'| > 2
                                                then 3
                                                else mb (3 + 1) z' (x, y)

So when we simplfy we find that the if branch is always the one chosen:

mb0 z (x, y) = let z' = z*z + x + y*i
               in  if |z'| > 2
                    then 0
                    else let z'' = z'*z' + x + y*i
                         in  if |z''| > 2
                             then 1
                             else let z''' = z''*z'' + x + y*i
                                  in  if |z'''| > 2
                                  then 2
                                  else 0

And our work here is done. Now we can throw away the definition of mb and always call mb0 instead. If we compile our unrolled variant with GHC and enable optimisation then the compiler will be able to optimise better and, going back to the real Fraskell code, we see a speed improvement of 20-25%! Hoorah!

We can also perform some simplifications on the code. one of the easiest that provides a significant speed improvement (due to (^) not being strict in its first argument) is transforming e^2 to let x = e in x * x. Between this and the unrolling, hand tweaked code can get a more than order of magnitude improvement.

However, things are not ideal. It takes us some effort to produce this unrolling and perform this simplification, and the resulting code is harder to read. Worse still, changing the bound on iterations means we now have to make drastic changes to our code. What would be nice is if we could get the speed of the unrolled version with the clarity of code of the original.

Template Haskell

Enter meta-programming! If you whisper the magic words (OK, OK, that's "-DTEMPLATE_HASKELL") while compiling Fraskell with a recent CVS snapshot of GHC then rather than the function being compiled as normal it gets intercepted by Template Haskell.

This provides us with the abstract syntax tree of the function, so we can then use normal Haskell functions to automate the unrolling before splicing the new definition back in. The actual transformation used isn't exactly as explained above as it proves to be easier to implement a slightly different transformation. For the implementation details see the source code and documentation, linked to below.

Some numbers

So how well does it work? These numbers, from generating a 1152x864 image on an AMD Athlon(tm) XP 2200+, show an almost order of magnitude speedup:

Building non-TH version

real    0m17.171s
user    0m16.440s
sys     0m0.720s


Running non-TH version

real    0m27.045s
user    0m26.890s
sys     0m0.110s


Building TH version

real    0m29.685s
user    0m28.310s
sys     0m1.240s


Running TH version

real    0m3.571s
user    0m3.420s
sys     0m0.130s


Comparing outputs


All done!

The code and documentation

Finally, the code (Fraskell 0.2.0) is available (released under the GNU GPL v2) if you want to have a play, as is the documentation of the Template Haskell utilities, including the Unroll and Simplify modules.

Should you still want it, the Fraskell 0.1.0 release is also still available (also released under the GNU GPL v2).


Back to the main page.

Random Image
Random Image
Random Image