Compiling a Mandelbrot program
Terry Robb
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 && kmandel, 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 && kmandel,
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
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
#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