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.
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.
|
