Compiling a Mandelbrot program

Terry Robb
<tdr@vaxc.cc.monash.edu.au>

Abstract

Here is a plot (with 128 by 128 resolution) of a particular part of the 
Mandelbrot set.

ListDensityPlot[ans//Transpose,
                ColorFunction->Hue, Mesh->False]; 1;



To produce this plot an array of integers called ans was first computed by 
various techniques. This Notebook looks at five different methods for 
computing ans.

The first method is to write a Function in Mathematica and then execute that 
function Ð this method is very slow and takes about half an hour to execute. 
The second and third method is to compile the function with either the 
Mathematica or InterCall compiler Ð this speeds up execution time by a factor 
of 20 and 60 respectively. The fourth and fifth method is to either Import 
with InterCall, or to Install with MathLink, a C-program Ð in either case the 
speed up in execution time is a factor of 180, i.e about 10 seconds. The 
timings for each of the methods are:

Mathematica non-compiled Function: 1800 Second of real-time
Mathematica Compiled Function:       92 Second of real-time
InterCall compiled Function:         34 Second of real-time
InterCall Imported C-program:        11 Second of real-time
MathLink Installed C-program:        10 Second of real-time

Compiling

The InterCall compiler

Consider a Mathematica function that can compute the Mandelbrot set within a 
radius r of some complex point c0 with the answer (that is the number of 
iterations, less than lim, until divergence) stored in an n by n array called 
ans. Here is such a function

mandel = Function[{ans,n,c0,r,lim},
  Module[{dx=r/(n-1),dy=r/(n-1),z=I,c=I,k=1},
    ans = Table[ z=0; c=c0+(2*i-1-n)*dx+(2*j-1-n)*I*dy; k=0;
                 While[Abs[z]<2 && k<lim, z=z^2+c; k++]; k,
                {i,1,n}, {j,1,n} ]], {HoldFirst}];

Calling the above function to compute the Mandelbrot set is quite slow 
however. For example

n=32; ans=.; mandel[ans, n, .37+.1*I, 0.005, 60];

will take about two minutes of real-time to execute. This time can be greatly 
decreased by using the compiler / interpreter available in InterCall version 
2.0.

To compile the mandel function we can use InterCall's AddDefault command. For 
this example the syntax would be

 AddDefault[$MANDEL->mandel, S[I[#2,#2],I,Z,R,I]];

where $MANDEL is just some dummy symbol that has been chosen to represent the 
compiled code. The last argument to AddDefault specifies the datatype of that 
code in a syntax described in the InterCall manual. This InterCall compilation 
takes about one second of cpu-time and the compiled code, called $MANDEL, is 
stored externally outside Mathematica where it would normally only be accessed 
by some other external program.

For our purpose we in fact want to access the code ourselves, and so we need 
to get some handle to $MANDEL that can be used at Mathematica's end. This 
handle can be obtained with InterCall's Peek command

mandel2 = Peek[$MANDEL];

After defining the above, we can call the InterCall compiled Mandelbrot code 
using the name mandel2. There is some subtlety about how mandel2 should be 
called however because the first argument to mandel2, internally refered to as 
ans, is a output-variable whose final value will be set by the external code. 
To pass such an object we just use the symbol None which indicates an empty 
initial value. Thus

n=32; mandel2[None, n, .37+.1*I, .005, 60]

will produce a 32 by 32 array centred at the complex point .37+.1i using a 
radius of .005. The array, corresponding to ans, will only be defined 
externally however, but we can use InterCall's Peek command to get it back 
into Mathematica. Since we have not given a proper external name to the array 
we can only refer to it by its position in the argument list. The array is the 
first argument to mandel2, and the command

ans=Peek[1];

will pull that first argument back to Mathematica. (It also possible to give a 
proper name to the array Ð see the InterCall manual for details.)

To call the InterCall compiled Mandelbrot code we thus type

n=32; mandel2[None, n, .37+.1*I, .005, 60]; ans=Peek[1];

and this takes just over two seconds to evaluate, compared with nearly two 
minutes for the equivalent non-compiled code.

The Mathematica compiler

We now compare the InterCall compiler / interpreter with the Mathematica 
compiler / interpreter. 

The Mathematica compiler does not handle vector (or array) arguments, and so 
for the Mandelbrot function we just compile the scalar operations and then use 
Table to evaluate that compiled code at each point. (There is some time 
penalty involved in doing this, but it is not enough to affect the conclusion 
that the InterCall interpreter is much faster than the Mathematica 
interpreter.) To compile the scalar part we use

mandel1 = Compile[{{n,_Integer},{c0,_Complex},{r,_Real},
                   {lim,_Integer},{i,_Integer},{j,_Integer}},
   Module[{dx=r/(n-1),dy=r/(n-1),z=I,c=I,k=1},
      z=0; c=c0+(2*i-1-n)*dx+(2*j-1-n)*I*dy; k=0;
      While[Abs[z]<2 && k<lim, z=z^2+c; k++]; k] ];

and then

n=32; ans=Table[mandel1[n, .37+.1*I, .005, 60, i, j],
                {i,1,n}, {j,1,n}];

can be used to evaluate the Mandelbrot set equivalent to the InterCall example 
used earlier.

The timing for the Mathematica compiled code above is about six seconds, which 
is about three times longer than the InterCall compiled code.

Importing vs Installing

Importing external code

InterCall can Import external code.

To do this you first use the AddDefault command to set up a default setting 
for the external code.

Once a default is set up, it can be inspected with the GetDefault command and 
modified with the SetDefault command. Re-issuing an AddDefault command will 
let you change the datatype associated with the routine Ð you don't need to re-
import to change the datatype.

Installing external code

MathLink can Install external code.

To do this you first write a MathLink template file to set up a link to the 
external code. You will also need to modify the C-source code and compile it 
with MathLink's mcc command.

Discussion

In the C-programs in this Notebook we have used hypot(zRE,zIM)<2 to test for 
divergence. It would be more efficient to test instead zRE*zRE+zIM*zIM<4 which 
would save doing an implicit square root, but for comparison with the 
Mathematica function mandel (which uses the test Abs[z]<2) we have not made 
that change.

The InterCall mandel3 C-program is virtually a literal translation of the 
Mathematica function mandel whereras the MathLink mandel4 C-program needs to 
be specially written with MathLink commands and a special template file. With 
MathLink programs, one needs to write at the C source level, whereas with 
InterCall all one needs is an existing object file or library archive.

InterCall's AddDefault command is the equivalent of MathLink's template file 
mechanism.

With InterCall one can dynamically change the AddDefault setting (with 
SetDefault), but with MathLink if one changes the template file then one has 
to recompile the C-program.

Various methods

This section contains the five different methods for computing the ans in the 
abstract.

mandel (Mathematica  Function)

Calling a Mathematica Function

Here is Mathematica Function called mandel

mandel = Function[{ans,n,c0,r,lim},
  Module[{dx=r/(n-1),dy=r/(n-1),z=I,c=I,k=1},
    ans = Table[ z=0; c=c0+(2*i-1-n)*dx+(2*j-1-n)*I*dy; k=0;
                 While[Abs[z]<2 && k<lim, z=z^2+c; k++]; k,
                {i,1,n}, {j,1,n} ]], {HoldFirst}];

To call the code type

n=128; ans=.; mandel[ans, n, .37+.1*I, .005, 60];

It takes about half an hour of real-time to call mandel as above. This time is 
quadratic with n.

mandel1 (Mathematica compiler/interpreter)

Compiling a Mathematica Function.

Here is a Mathematica Compiled Function called mandel1

mandel1 = Compile[{{n,_Integer},{c0,_Complex},{r,_Real},
                   {lim,_Integer},{i,_Integer},{j,_Integer}},
   Module[{dx=r/(n-1),dy=r/(n-1),z=I,c=I,k=1},
      z=0; c=c0+(2*i-1-n)*dx+(2*j-1-n)*I*dy; k=0;
      While[Abs[z]<2 && k<lim, z=z^2+c; k++]; k] ];

To call the code type

n=128; ans=Table[mandel1[n, .37+.1*I, .005, 60, i, j],
                 {i,1,n}, {j,1,n}]; 1;

It takes over 90 seconds of real-time to call mandel1 as above.

Here is the Mathematica pseudo-code associated with mandel1

mandel1//InputForm

CompiledFunction[{_Integer, _Complex, _Real, _Integer, 
   _Integer, _Integer}, {3, 27, 17, 12}, 
  {{1, 17}, {3, 1, 0}, {5, 2, 0}, {4, 3, 0}, {3, 4, 1}, 
   {3, 5, 2}, {3, 6, 3}, {12, -1, 4}, {32, 0, 4, 5}, 
   {20, 5, 1}, {41, 1, 2}, {36, 0, 2, 3}, {12, -1, 6}, 
   {32, 0, 6, 7}, {20, 7, 4}, {41, 4, 5}, {36, 0, 5, 6}, 
   {14, 0. + 1.*I, 1}, {14, 0. + 1.*I, 2}, {12, 1, 8}, 
   {12, 0, 9}, {20, 9, 7}, {13, 0., 8}, {21, 7, 8, 3}, 
   {18, 3, 1}, {12, 2, 10}, {35, 10, 2, 11}, 
   {12, -1, 12}, {35, 0, 13}, {38, 13, 14}, 
   {32, 11, 12, 14, 15}, {20, 15, 9}, {36, 9, 3, 10}, 
   {12, 2, 16}, {35, 16, 3, 17}, {12, -1, 18}, 
   {35, 0, 19}, {38, 19, 20}, {32, 17, 18, 20, 21}, 
   {14, 0. + 1.*I, 4}, {20, 21, 11}, {13, 0., 12}, 
   {21, 11, 12, 5}, {13, 0., 13}, {21, 6, 13, 6}, 
   {37, 5, 4, 6, 7}, {13, 0., 14}, {21, 10, 14, 8}, 
   {34, 0, 8, 7, 9}, {18, 9, 2}, {12, 0, 22}, 
   {16, 22, 8}, {67, 1, 15}, {12, 2, 23}, {20, 23, 16}, 
   {75, 15, 16, 0}, {74, 8, 1, 1}, {78, 0, 1, 2}, 
   {69, 2, 9}, {37, 1, 1, 10}, {34, 10, 2, 11}, 
   {18, 11, 1}, {16, 8, 24}, {12, 1, 25}, 
   {32, 8, 25, 26}, {16, 26, 8}, {70, -14}, {7, 8}}, 
  Function[{n, c0, r, lim, i, j}, 
   Module[{dx = r/(n - 1), dy = r/(n - 1), z = I, c = I, 
     k = 1}, z = 0; c = 
      c0 + (2*i - 1 - n)*dx + (2*j - 1 - n)*I*dy; k = 0; 
     While[Abs[z] < 2 && k < lim, z = z^2 + c; k++]; k]]]

mandel2 (InterCall compiler/interpreter)

Compiling a Function using InterCall.

This example needs InterCall switched on

Needs["InterCall`"];
InterCall[On];

Here is an InterCall compiled Function called mandel2

mandel = Function[{ans,n,c0,r,lim},
  Module[{dx=r/(n-1),dy=r/(n-1),z=I,c=I,k=1},
    ans = Table[ z=0; c=c0+(2*i-1-n)*dx+(2*j-1-n)*I*dy; k=0;
                 While[Abs[z]<2 && k<lim, z=z^2+c; k++]; k,
                {i,1,n}, {j,1,n} ]], {HoldFirst}];
mandel2 = Peek @ AddDefault[$MANDEL->mandel,
                            S[I[#2,#2],I,Z,R,I]];

To call the code type

n=128; mandel2[None, n, .37+.1*I, .005, 60]; ans=Peek[1]; 1;

It takes about 30 seconds of real-time to call mandel2 as above, and then 
about 3 seconds of real-time to Peek back the ans. These times are both 
quadratic with n.

Here is the InterCall pseudo-code associated with mandel2

mandel2//InputForm

ExternalCode[1, S[I[$Arg[1, 2], $Arg[1, 2]], I, Z, R, I], 
  1, RemoteFunction[{_Integer, _Integer, _Complex, _Real, 
    _Integer}, {5, 63, 25, 12}, 
   {{1, 17}, {3, 1, 0}, {3, 2, 1}, {5, 3, 0}, {4, 4, 0}, 
    {3, 5, 2}, {12, -1, 3}, {32, 1, 3, 4}, {20, 4, 1}, 
    {41, 1, 2}, {36, 0, 2, 3}, {12, -1, 5}, 
    {32, 1, 5, 6}, {20, 6, 4}, {41, 4, 5}, {36, 0, 5, 6}, 
    {14, 0. + 1.*I, 1}, {14, 0. + 1.*I, 2}, {12, 1, 7}, 
    {12, -1, 8}, {16, 8, 9}, {12, 1, 10}, {16, 10, 11}, 
    {12, 1, 12}, {16, 12, 13}, {38, 11, 14}, 
    {32, 1, 14, 15}, {20, 13, 7}, {41, 7, 8}, 
    {20, 15, 9}, {36, 9, 8, 10}, {63, 10, 16}, 
    {16, 16, 17}, {16, 11, 18}, {12, 1, 19}, 
    {32, 9, 19, 20}, {16, 20, 9}, {35, 9, 13, 21}, 
    {32, 11, 21, 22}, {16, 22, 18}, {76, 9, 17, 0}, 
    {69, 0, 81}, {12, -1, 23}, {16, 23, 24}, {12, 1, 25}, 
    {16, 25, 26}, {12, 1, 27}, {16, 27, 28}, 
    {38, 26, 29}, {32, 1, 29, 30}, {20, 28, 11}, 
    {41, 11, 12}, {20, 30, 13}, {36, 13, 12, 14}, 
    {63, 14, 31}, {16, 31, 32}, {16, 26, 33}, 
    {12, 1, 34}, {32, 24, 34, 35}, {16, 35, 24}, 
    {35, 24, 28, 36}, {32, 26, 36, 37}, {16, 37, 33}, 
    {76, 24, 32, 1}, {69, 1, -30}, {12, 0, 38}, 
    {12, 0, 39}, {20, 39, 15}, {13, 0., 16}, 
    {21, 15, 16, 3}, {18, 3, 1}, {12, 2, 40}, 
    {35, 40, 18, 41}, {12, -1, 42}, {35, 1, 43}, 
    {38, 43, 44}, {32, 41, 42, 44, 45}, {20, 45, 17}, 
    {36, 17, 3, 18}, {12, 2, 46}, {35, 46, 33, 47}, 
    {12, -1, 48}, {35, 1, 49}, {38, 49, 50}, 
    {32, 47, 48, 50, 51}, {14, 0. + 1.*I, 4}, 
    {20, 51, 19}, {13, 0., 20}, {21, 19, 20, 5}, 
    {13, 0., 21}, {21, 6, 21, 6}, {37, 5, 4, 6, 7}, 
    {13, 0., 22}, {21, 18, 22, 8}, {34, 0, 8, 7, 9}, 
    {18, 9, 2}, {12, 0, 52}, {16, 52, 7}, {67, 1, 23}, 
    {12, 2, 53}, {20, 53, 24}, {75, 23, 24, 2}, 
    {74, 7, 2, 3}, {78, 2, 3, 4}, {69, 4, 9}, 
    {37, 1, 1, 10}, {34, 10, 2, 11}, {18, 11, 1}, 
    {16, 7, 54}, {12, 1, 55}, {32, 7, 55, 56}, 
    {16, 56, 7}, {70, -14}, {16, 7, 38}, {12, -1, 57}, 
    {32, 33, 57, 58}, {35, 1, 58, 59}, {32, 18, 59, 60}, 
    {0}, {83, 0, 60, 38}, {16, 38, 62}, {70, -64}, {10}}, 
   Function[{ans, n, c0, r, lim}, 
    Module[{dx = r/(n - 1), dy = r/(n - 1), z = I, c = I, 
      k = 1}, ans = 
      Table[z = 0; c = 
         c0 + (2*i - 1 - n)*dx + (2*j - 1 - n)*I*dy; 
        k = 0; While[Abs[z] < 2 && k < lim, 
         z = z^2 + c; k++]; k, {i, 1, n}, {j, 1, n}]]]]]

mandel3 (InterCall C-program)

Importing a program using InterCall.

This example needs InterCall

Needs["InterCall`"];

Here is the C source code associated with mandel3

!!mandel3.c

/* % cc -O -c mandel3.c */
#include <math.h>

void mandel3_(ANS,N,C0,R,LIM) long *ANS; long *N; double *C0, *R; long *LIM;
{
  double dx=(*R)/((*N)-1), dy=(*R)/((*N)-1), zRE,zIM, cRE,cIM, tmp;
  long i,j,k;
  for (i=1; i<=(*N); i++) {
    for (j=1; j<=(*N); j++) {
      zRE = 0; cRE = C0[0] + (2*i-1-(*N))*dx;
      zIM = 0; cIM = C0[1] + (2*j-1-(*N))*dy;
      k   = 0;
      while (hypot(zRE,zIM)<2 && k<(*LIM)) {
        tmp = 2*zRE*zIM + cIM;
        zRE = zRE*zRE - zIM*zIM + cRE;
        zIM = tmp;
        k++; }
      ANS[(i-1)+(j-1)*(*N)] = k;
    }
  }
}

To import mandel3 type

AddDefault[ MANDEL3[$ANS:>Out,$N,$C0,$R,$LIM],
            S[I[$N,$N],I,Z,R,I],
            "mandel3.o" ];
Import[{MANDEL3}];

To call the code type

n=128; ans = mandel3[n, .37+.10*I, .005, 60]; 1;

It takes about 11 seconds of real-time to call mandel3 as above.

mandel4 (MathLink C-program)

Installing a program using MathLink.

Here is the C source code associated with mandel4

!!mandel4.tm

/* % mcc -O -o mandel4 mandel4.tm */
:Begin:
:Function:      mandel4
:Pattern:       mandel4[n_Integer,c0_Complex,r_Real,lim_Integer]
:Arguments:     {n,Re[c0],Im[c0],r,lim}
:ArgumentTypes: {Integer,Real,Real,Real,Integer}
:ReturnType:    Manual
:End:

#include <math.h>
#include "mathlink.h"

int main(argc,argv) int argc; char* argv[]; { return MLMain(argc,argv); }

void mandel4(n,c0RE,c0IM,r,lim) long n; double c0RE,c0IM,r; long lim;
{
  double dx=r/(n-1), dy=r/(n-1), zRE,zIM, cRE,cIM, tmp;
  long i,j,k;

  MLPutFunction(stdlink,"List",n);
  for (i=1; i<=n; i++) {
    MLPutFunction(stdlink,"List",n);
    for (j=1; j<=n; j++) {
      zRE = 0; cRE = c0RE + (2*i-1-n)*dx;
      zIM = 0; cIM = c0IM + (2*j-1-n)*dy;
      k   = 0;
      while (hypot(zRE,zIM)<2 && k<lim) {
        tmp = 2*zRE*zIM + cIM;
        zRE = zRE*zRE - zIM*zIM + cRE;
        zIM = tmp;
        k++; }
      MLPutInteger(stdlink,k);
    }
  }
}

To install mandel4 type

Install["mandel4"];

To call the code type

n=128; ans = mandel4[n, .37+.1*I, .005, 60]; 1;

It takes about 10 seconds of real-time to call mandel4 as above.

Timings

The interpreter associated with the InterCall compiler is faster than the 
interpreter associated with the Mathematica compiler. For the Mandelbrot 
example the InterCall interpreter is about three times faster, although for 
most functions it is typically only about twice as fast. Here are the timings 
for the five different methods of solving a Mandelbrot problem on a 128 by 128 
grid.

             n=128; c0=.37+.1I; r=.005; lim=60;
Mathematica non-compiled Function: 1800   Second of real-time
Mathematica Compiled Function:       92.3 Second of real-time
InterCall compiled Function:         33.8 Second of real-time
InterCall Imported C-program:        10.9 Second of real-time
MathLink Installed C-program:         9.9 Second of real-time

The MathLink C-program is slightly faster than the InterCall C-program because 
InterCall has a higher initial overhead when calling external routines due to 
it having to make decisions about default actions, and also due to the manner 
in which communication has to be done. Commands and timings for a Mandelbrot 
problem on a 32 by 32 grid are given below.

Mathematica non-compiled Function

n=32; ans=.; mandel[ans, n, .37+.1*I, .005, 60];
     (* ShowTime::timing: 98.6 Second (111.04352) *)

Mathematica Compiled Function

n=32; ans=Table[mandel1[n, .37+.1*I, .005, 60, i, j],
                {i,1,n}, {j,1,n}];
     (* ShowTime::timing: 5.08333 Second (5.41011) *)

InterCall compiled Function

n=32; mandel2[None, n, .37+.1*I, .005, 60]; ans=Peek[1];
     (* ShowTime::timing: 0.233333 Second (2.12875) *)

InterCall Imported C-program

n=32; ans=mandel3[n, .37+.1*I, .005, 60];
     (* ShowTime::timing: 0.2 Second (0.81642) *)

MathLink Installed C-program

n=32; ans=mandel4[n, .37+.1*I, .005, 60];
     (* ShowTime::timing: 0.0833333 Second (0.63277) *)

Conclusion

Calling InterCall imported or MathLink installed C-code is about three times 
faster than calling InterCall compiled code, and calling InterCall compiled 
code is about three times faster than calling Mathematica compiled code.

Apart from speed the InterCall compiler has various advantages over its 
Mathematica counterpart. 

Also importing a C-program with InterCall has advantages over installing with 
MathLink.

References

InterCall

InterCall was developed by Terry Robb. For more information on InterCall send 
email to intercall_forward@wri.com. An order form is available via MathSource Ð
 email the line send 0202-587-0011 to mathsource@wri.com for details. 
Information can also be posted or faxed to Analytica below.

Analytica                     Fax: +61 (9) 386 5666
P.O. Box 522
Nedlands 6009
Perth WA
AUSTRALIA.