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:

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

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

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.

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.