28
$\begingroup$

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:

enter image description here

Improve this question
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – user4324
    Commented Oct 19, 2012 at 3:42
  • $\begingroup$ 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! $\endgroup$
    – Verbeia
    Commented Oct 19, 2012 at 5:55

5 Answers 5

Reset to default
27
$\begingroup$

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}]]

enter image description here

Here's a smoother version with 0.5 as the nearness limit:

smoother

And some variations:

variation1 variation2

And some animations of how the image gets constructed:

Individual gaussian integers Sum over gaussian integers

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.

Improve this answer
$\endgroup$
6
  • 1
    $\begingroup$ Nice image, maybe the final ] fell off the Image[...] during copying? $\endgroup$
    – image_doctor
    Commented Oct 18, 2012 at 10:59
  • 3
    $\begingroup$ @image_doctor I've added some more images, and 8.3% more right brackets! $\endgroup$
    – wxffles
    Commented Oct 18, 2012 at 19:20
  • $\begingroup$ Were the last few ones done with ReliefPlot[]? $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Oct 18, 2012 at 22:34
  • $\begingroup$ 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. $\endgroup$
    – wxffles
    Commented Oct 18, 2012 at 22:43
  • 1
    $\begingroup$ Finally the question becomes interesting! But what precisely is being varied in the animations? $\endgroup$
    – whuber
    Commented Oct 19, 2012 at 14:12
13
$\begingroup$

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:)

Example fractal

Improve this answer
$\endgroup$
12
$\begingroup$

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]]

Pseudo fractal

Improve this answer
$\endgroup$
2
  • 3
    $\begingroup$ 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]}]]] $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Oct 17, 2012 at 16:31
  • $\begingroup$ @J.M. Nice! That increases the speed by 50%. $\endgroup$
    – whuber
    Commented Oct 17, 2012 at 16:34
6
$\begingroup$

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.

Improve this answer
$\endgroup$
5
  • $\begingroup$ a != 0 || b != 0 could be changed into ! (a == 0 && b == 0). $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Oct 17, 2012 at 12:50
  • $\begingroup$ @J.M. I don't see much difference. Should I? $\endgroup$
    – Mr.Wizard
    Commented Oct 17, 2012 at 12:52
  • $\begingroup$ 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; $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Oct 17, 2012 at 12:55
  • $\begingroup$ @J.M. or leave it out entirely and just use Quiet, which seems to work as well. $\endgroup$
    – Mr.Wizard
    Commented Oct 17, 2012 at 12:57
  • $\begingroup$ @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. $\endgroup$
    – jVincent
    Commented Oct 17, 2012 at 13:58
6
$\begingroup$

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]

plot of check array

Why it looks nothing like the picture in the OP, I don't know...

Improve this answer
$\endgroup$
3
  • $\begingroup$ Clean, but slower. $\endgroup$
    – Mr.Wizard
    Commented Oct 17, 2012 at 12:55
  • $\begingroup$ Unfortunately, I can't get ParallelTable[] to work on my end, so someone might want to see what happens if Table[] is replaced accordingly. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Oct 17, 2012 at 12:57
  • $\begingroup$ I found Table just as fast as ParallelTable with default options but I didn't try different Methods. $\endgroup$
    – Mr.Wizard
    Commented Oct 17, 2012 at 12:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.