# Making fractals with Mathematica

I recently saw this post on math.stackexchange and was curious as to how to generate the image in Mathematica. I tried the following naive approach; however, it is extremely slow.

Clear[check, GaussianIntegerQ]
GaussianIntegerQ[a_] := If[IntegerQ[Re[a]] && IntegerQ[Im[a]], True, False]
check[a_] := Block[{d = 0},Do[If[GaussianIntegerQ[c (1 + I)/a], d++], {c, 1, 100}]; d];
ArrayPlot[ParallelTable[If[a != 0 || b != 0, check[a + b I], 0], {a, -1, 1, 1/100},
{b,-1, 1, 1/100}], ColorFunction -> (GrayLevel[#] &)] // AbsoluteTiming

(*{22.5931794, img}*)


I tried making it faster, but the speed-up wasn't much:

Clear[check]
check[a_, b_] := Block[{d = 0},Do[If[IntegerQ[(a c + b c)/(a^2 + b^2)] &&
IntegerQ[(a c - b c)/(a^2 + b^2)], d++], {c, 1, 100}]; d]
ArrayPlot[ParallelTable[If[a != 0 || b != 0, check[a, b], 0], {a, -1, 1,1/100},
{b, -1, 1, 1/100}], ColorFunction -> (GrayLevel[#] &)] // AbsoluteTiming

(*{15.5660219, img}*)


Could anyone offer suggestions on how to make it faster? (for what it's worth, here is C-code from a comment on the blog post)

The final result should look something like: • In addition to the compilation people have suggested, you can take advantage of the symmetry and only generate one quadrant of it. Then flip that quadrant around to display the whole image. On my laptop that was faster than compilation alone.
– user4324
Oct 19, 2012 at 3:42
• Dear @RichardTodd, welcome to Mathematica.SE! You have made a good suggestion, but it wasn't really enough to be a separate answer. I know this might have been because you don't have enough reputation to comment yet. If you have code you can post, to complement your comment, you could post that as an answer, together with the explanatory text that I have converted into a comment. Once again, welcome and thanks for your contribution! Oct 19, 2012 at 5:55

The whole "fractal" is an exercise in rounding errors. Following all the links to some code, we find that something is considered an integer if its fractional part is less than 0.1. Using something similar to Mr.Wizard's answer:

inQ = Abs[FractionalPart[N[#, 16]]] < 0.1 &;
check[0 | 0., 0 | 0.] := 0;
check[a_, b_] :=
With[{p = (a + b)/(a^2 + b^2), q = (a - b)/(a^2 + b^2)},
Sum[Boole[inQ[c p] && inQ[c q]], {c, 100}]];

Image[Table[(0.01 #)^(1/4) &@ check[a, b], {a, -1, 1, 0.0025}, {b, -1, 1, 0.0025}]] Here's a smoother version with 0.5 as the nearness limit: And some variations:  And some animations of how the image gets constructed:  Edit for the curious:

The left animates the binary images you get from considering each gaussian integer individually: $1+i$, then $2+2i$, etc. The images on the right are the sums from $1+i$ to $k+ki$, or essentially the sums of the binary images on the left. Also the range is -5 to 5 instead of the -1 to 1 of the original.

• Nice image, maybe the final ] fell off the Image[...] during copying? Oct 18, 2012 at 10:59
• @image_doctor I've added some more images, and 8.3% more right brackets! Oct 18, 2012 at 19:20
• Were the last few ones done with ReliefPlot[]? Oct 18, 2012 at 22:34
• They were done with Abs[FractionalPart[a]] < 0.5], at high resolution, then gamma corrected and scaled down. I highly recommend Ajasja's compiled version unless you have a lot of time on your hands. Oct 18, 2012 at 22:43
• Finally the question becomes interesting! But what precisely is being varied in the animations? Oct 19, 2012 at 14:12

This is a compiled version of @wxffles answer (since I got bored waiting for the uncompiled version to finish on my slow home computer:)

inQ = Compile[{a}, Abs[FractionalPart[a]] < 0.1];
check = Compile[{a, b},
With[{p = (a + b)/(a^2 + b^2), q = (a - b)/(a^2 + b^2)},
Sum[Boole[inQ[c p] && inQ[c q]], {c, 100}]],
CompilationOptions -> {"InlineCompiledFunctions" -> True,
"InlineExternalDefinitions" -> True}];

getFractal =
Compile[{},
Table[(0.01 #)^(1/4) &@check[a, b], {a, -1, 1, 0.0025}, {b, -1, 1,
0.0025}],
CompilationOptions -> {"InlineCompiledFunctions" -> True,
"InlineExternalDefinitions" -> True,
"ExpressionOptimization" -> True}, RuntimeOptions -> {"Speed"}];

AbsoluteTiming[Image[getFractal[]]]


The timings on my work computer (AMD Phenom II X6 1090T) are:

 Uncompiled           18.0 s
Compiled to to WVM   10.6 s
Compiled to C         1.4 s


So we get an order of magnitude speedup using compile. Oh, here is an image generated, so we see it is the same:) The problem with the speed is that check is using brute force to count the Gaussian integers among the first 100 multiples of $(1+i)/(a + bi)$ by enumerating and checking them all. This count can be computed directly for about two orders of magnitude speedup simply by finding the least common denominator of the real and imaginary parts of the quotient:

check[a_] := With[{u = (1 + I)/a, n = 100}, Floor[n / LCM @@ (Denominator /@ {Re[u], Im[u]})]];


Generating the image in the question takes 0.233 seconds on my machine and is identical to the original image (9.34 seconds), for about a forty fold speedup. I suspect Mathematica wizards could improve my crude implementation and perhaps wring another factor of two out of its performance.

With[{g = 1/10, n = 500}, ArrayPlot[t = ParallelTable[
If[a != 0 || b != 0, check[a + b I ], 0], {a, -1, 1, 1/n}, {b, -1, 1, 1/n}],
ColorFunction -> (GrayLevel[#^g] &), ImageSize -> n]] • Since Denominator[] is Listable, and Quotient[m, n] == Floor[m/n], one can have check[a_] := With[{u = (1 + I)/a, n = 100}, Quotient[n, LCM @@ Denominator[{Re[u], Im[u]}]]] Oct 17, 2012 at 16:31
• @J.M. Nice! That increases the speed by 50%. Oct 17, 2012 at 16:34

As a first pass this is about 40% faster on my machine:

iQ = # == Round@# &;

check[a_, b_] :=
With[{s = (a^2 + b^2)},
Sum[Boole[iQ[(a c + b c)/s] && iQ[(a c - b c)/s]], {c, 100}]]

ArrayPlot[
Table[If[a != 0 || b != 0, check[a, b], 0], {a, -1, 1, 1/100}, {b, -1, 1, 1/100}],
ColorFunction -> GrayLevel]


This is using machine-precision arithmetic.

• a != 0 || b != 0 could be changed into ! (a == 0 && b == 0). Oct 17, 2012 at 12:50
• @J.M. I don't see much difference. Should I? Oct 17, 2012 at 12:52
• It's your call; actually, now that I think about it, maybe it might be more expedient to add a definition for check[]: check[a : 0 | 0., b : 0 | 0.] = 0; Oct 17, 2012 at 12:55
• @J.M. or leave it out entirely and just use Quiet, which seems to work as well. Oct 17, 2012 at 12:57
• @Mr.Wizard It seems that using Compile this can be speed up even more dramatically, going from 14 seconds to 0.921 on my system. Oct 17, 2012 at 13:58

Try this:

SetAttributes[check, Listable]
check = 0;
check[a_] := Count[Divisible[(1 + I) Range, a], True];

ArrayPlot[check[Table[a + b I, {a, -1, 1, 1/100}, {b, -1, 1, 1/100}]],
ColorFunction -> GrayLevel] Why it looks nothing like the picture in the OP, I don't know...

• Clean, but slower. Oct 17, 2012 at 12:55
• Unfortunately, I can't get ParallelTable[] to work on my end, so someone might want to see what happens if Table[] is replaced accordingly. Oct 17, 2012 at 12:57
• I found Table just as fast as ParallelTable` with default options but I didn't try different Methods. Oct 17, 2012 at 12:59