Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

share|improve this question
    
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

5 Answers 5

up vote 21 down vote accepted

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.

share|improve this answer
1  
Nice image, maybe the final ] fell off the Image[...] during copying? –  image_doctor Oct 18 '12 at 10:59
1  
@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[]? –  J. M. 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
1  
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:)

Example fractal

share|improve this answer

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

share|improve this answer
3  
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]}]]] –  J. M. 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.

share|improve this answer
    
a != 0 || b != 0 could be changed into ! (a == 0 && b == 0). –  J. M. 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; –  J. M. 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]

plot of check array

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

share|improve this answer
    
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. –  J. M. 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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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