# 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 '12 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! –  Verbeia Oct 19 '12 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? –  image_doctor Oct 18 '12 at 10:59
@image_doctor I've added some more images, and 8.3% more right brackets! –  wxffles Oct 18 '12 at 19:20
Were the last few ones done with ReliefPlot[]? –  Guess who it is. Oct 18 '12 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. –  wxffles Oct 18 '12 at 22:43
Finally the question becomes interesting! But what precisely is being varied in the animations? –  whuber Oct 19 '12 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]}]]] –  Guess who it is. Oct 17 '12 at 16:31
@J.M. Nice! That increases the speed by 50%. –  whuber Oct 17 '12 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). –  Guess who it is. Oct 17 '12 at 12:50
@J.M. I don't see much difference. Should I? –  Mr.Wizard Oct 17 '12 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; –  Guess who it is. Oct 17 '12 at 12:55
@J.M. or leave it out entirely and just use Quiet, which seems to work as well. –  Mr.Wizard Oct 17 '12 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. –  jVincent Oct 17 '12 at 13:58

Try this:

SetAttributes[check, Listable]
check[0] = 0;
check[a_] := Count[Divisible[(1 + I) Range[100], 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. –  Mr.Wizard Oct 17 '12 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. –  Guess who it is. Oct 17 '12 at 12:57
I found Table just as fast as ParallelTable` with default options but I didn't try different Methods. –  Mr.Wizard Oct 17 '12 at 12:59