# Speediest Julia Set "By Hand"

While being aware of the built-in command, I would like to show my students how the Julia set is actually generated. The code I came up with is slow, even though I am quite happy with the picture it ultimately creates.

Can you help me with recommending code replacements that would make the image generation speedier?

\[Delta]=0.003;
c=-0.04505+0.66142I;
p=#^2+c&;
grayLevels=200;
iter[z_]:=Length[NestWhileList[p,z,Abs[#]<4&,1,grayLevels]]
range=Flatten[Table[a+b I,{a,-1.5,1.5,\[Delta]},{b,-1.3,1.3,\[Delta]}]];
results=Map[{{Re[#],Im[#]},iter[#]}&,range];
cRange=Union[results[[All,2]]];
cMax=Max[cRange];
cMin=Min[cRange];
cGray=Map[GrayLevel[(#-cMin)/(cMax-cMin)]&,cRange];
ListPlot[Table[Cases[results,x_/;x[[2]]==c][[All,1]],{c,cRange}],AspectRatio->Automatic,PlotStyle->cGray,Axes->False,ImageSize->Large,PlotMarkers->{Automatic, Tiny}]


• See here
– I.M.
Commented Mar 3, 2022 at 3:46
• Ah, thanks, I.M. One issue that seems to hamper the speed of my code is the rendering. In the link you provided, they only graph the result of the backwards iteration of the polynomial. I, however, test each point in the box for how fast it drifts off to infinity, so I need to plot each point with its appropriate color value. Indeed, if I separate "computation" and "rendering" into separate cells, it's the rendering that eats most of the time. I am baffled how the built-in command manages to render so quickly. Commented Mar 4, 2022 at 15:54

Plotting a million points is rather slow, instead you can use Image. This requires the results in a 2D array, so remove the Flatten from the definition of range. Something like this:

range = Table[a + b I, {a, -1.5, 1.5, δ}, {b, -1.3, 1.3, δ}];
results = Map[iter, range, {2}];
Image[results/grayLevels]


That makes the display fast, but the calculation is still slow. Compile is a good option here. To keep things simple I've hard coded the values of c and grayLevels into the compiled function, but you could easily make them input parameters.

This should run in a couple of seconds.

iter = Compile[{{z, _Complex}},
Block[{n = 0, x = z, c = -0.04505 + 0.66142 I},
While[n < 200 && Abs[x] < 4, x = x^2 + c; n++];
n/200], RuntimeAttributes -> {Listable}];

δ = 0.003;

range = Table[Complex[a, b], {a, -1.5, 1.5, δ}, {b, -1.3, 1.3, δ}];

Image[iter[range]]


If you want it really fast the new LLVM compiler does a good job. This is quick enough to put in a Manipulate and vary the c parameter in real time.

cf = FunctionCompile[
FunctionDeclaration[escapetime,
Typed[{"Complex128", "Complex128"} -> "Integer64"]@Function[{z, c},
Block[{n = 0, x = z},
While[n < 200 && Abs[x] < 4, x = x^2 + c; n++];
n]]],
Function[{Typed[n, "Integer64"], Typed[c, "Complex128"]},
Table[escapetime[Complex[a, b], c]/200.0,
{a, -1.5, 1.5, 3.0/n}, {b, -1.5, 1.5, 3.0/n}]]]

Manipulate[
Image[cf[250, Complex @@ pt]],
{{pt, {-0.04505, 0.66142}}, {-1, -1}, {1, 1}}]


• Darn, can't upvote a second time. Commented Mar 6, 2022 at 22:27
• The "image" one is still high level enough to be easily understandable without rigorous study of the manual, but it runs in about 0.6 seconds on my computer. It's PERFECT! Commented Mar 8, 2022 at 14:55
• +1 for the FunctionCompile showcase! However on my computer (win10+M13.0), cf performs slightly worse than iter when I set options Parallelization->True, RuntimeAttributes->{Listable}, RuntimeOptions->"Speed", CompilationTarget->"C". Does FunctionCompile in current version really have performance advantage over the old Compile in certain cases? Commented Jun 23, 2022 at 14:53
• @Silvia, I think I would expect broadly similar performance, at least for simple calculations like this. In my view the main advantage of FunctionCompile right now is that it works out of the box without needing to install a separate compiler. Commented Jun 25, 2022 at 22:33
• ParallelParallelTable instead of Table inside FunctionCompile` seems to give some extra boost. Commented Jul 6, 2022 at 15:04