Given three mutually tangent circles, we can always find two more circles which are tangent to all three of those. These two are called Apollonian circles. Note that one of the Apollonian circles might actually be around the three initial circles.

Starting from three tangent circles, we can create a fractal called an Apollonian gasket, by the following process:

1. Call the initial 3 circles the parent circles
2. Find the parent circles' two Apollonian circles
3. For each Apollonian circle:
1. For each pair of the three pairs of parent circles:
1. Call the Apollonian circle and the two parent circles the new set of parent circles and start over from step 2.

E.g. starting with circles of equal size, we get:

Image found on Wikipedia

There's one more bit of notation we need. If we have a circle of radius r with centre (x, y), we can define it's curvature as k = ±1/r. Usually k will be positive, but we can use negative k to denote the circle that encloses all the other circles in the gasket (i.e. all tangents touch that circle from the inside). Then we can specify a circle with a triplet of numbers: (k, x*k, y*k).

For the purpose of this question, we will assume positive integer k and rational x and y.

Further examples for such circles can be found in the Wikipedia article.

The Challenge

You will be given 4 circle specifications, each of which will look like (14, 28/35, -112/105). You can use any list format and division operator that is convenient, such that you can simply eval the input if you wish to. You may assume that the 4 circles are indeed tangent to each other, and that the first of them has negative curvature. That means you are already given the surrounding Apollonian circle of the other three. For a list of valid example inputs, see the bottom of the challenge.

Write a program or function which, given this input, draws an Apollonian gasket.

You may take input via function argument, ARGV or STDIN and either render the fractal on screen or write it to an image file in a format of your choice.

If the resulting image is rasterised, it must be at least 400 pixels on each side, with less than 20% padding around the largest circle. You may stop recursing when you reach circles whose radius is less than a 400th of the largest input circle, or circles which are smaller than a pixel, whichever happens first.

You must draw only circle outlines, not full discs, but the colours of background and lines are your choice. The outlines must not be wider than a 200th of the outer circles diameter.

This is code golf, so the shortest answer (in bytes) wins.

Example Inputs

Here are all integral gaskets from the Wikipedia article converted to the prescribed input format:

[[-1, 0, 0], [2, 1, 0], [2, -1, 0], [3, 0, 2]]
[[-2, 0, 0], [3, 1/2, 0], [6, -2, 0], [7, -3/2, 2]]
[[-3, 0, 0], [4, 1/3, 0], [12, -3, 0], [13, -8/3, 2]]
[[-3, 0, 0], [5, 2/3, 0], [8, -4/3, -1], [8, -4/3, 1]]
[[-4, 0, 0], [5, 1/4, 0], [20, -4, 0], [21, -15/4, 2]]
[[-4, 0, 0], [8, 1, 0], [9, -3/4, -1], [9, -3/4, 1]]
[[-5, 0, 0], [6, 1/5, 0], [30, -5, 0], [31, -24/5, 2]]
[[-5, 0, 0], [7, 2/5, 0], [18, -12/5, -1], [18, -12/5, 1]]
[[-6, 0, 0], [7, 1/6, 0], [42, -6, 0], [43, -35/6, 2]]
[[-6, 0, 0], [10, 2/3, 0], [15, -3/2, 0], [19, -5/6, 2]]
[[-6, 0, 0], [11, 5/6, 0], [14, -16/15, -4/5], [15, -9/10, 6/5]]
[[-7, 0, 0], [8, 1/7, 0], [56, -7, 0], [57, -48/7, 2]]
[[-7, 0, 0], [9, 2/7, 0], [32, -24/7, -1], [32, -24/7, 1]]
[[-7, 0, 0], [12, 5/7, 0], [17, -48/35, -2/5], [20, -33/35, 8/5]]
[[-8, 0, 0], [9, 1/8, 0], [72, -8, 0], [73, -63/8, 2]]
[[-8, 0, 0], [12, 1/2, 0], [25, -15/8, -1], [25, -15/8, 1]]
[[-8, 0, 0], [13, 5/8, 0], [21, -63/40, -2/5], [24, -6/5, 8/5]]
[[-9, 0, 0], [10, 1/9, 0], [90, -9, 0], [91, -80/9, 2]]
[[-9, 0, 0], [11, 2/9, 0], [50, -40/9, -1], [50, -40/9, 1]]
[[-9, 0, 0], [14, 5/9, 0], [26, -77/45, -4/5], [27, -8/5, 6/5]]
[[-9, 0, 0], [18, 1, 0], [19, -8/9, -2/3], [22, -5/9, 4/3]]
[[-10, 0, 0], [11, 1/10, 0], [110, -10, 0], [111, -99/10, 2]]
[[-10, 0, 0], [14, 2/5, 0], [35, -5/2, 0], [39, -21/10, 2]]
[[-10, 0, 0], [18, 4/5, 0], [23, -6/5, -1/2], [27, -4/5, 3/2]]
[[-11, 0, 0], [12, 1/11, 0], [132, -11, 0], [133, -120/11, 2]]
[[-11, 0, 0], [13, 2/11, 0], [72, -60/11, -1], [72, -60/11, 1]]
[[-11, 0, 0], [16, 5/11, 0], [36, -117/55, -4/5], [37, -112/55, 6/5]]
[[-11, 0, 0], [21, 10/11, 0], [24, -56/55, -3/5], [28, -36/55, 7/5]]
[[-12, 0, 0], [13, 1/12, 0], [156, -12, 0], [157, -143/12, 2]]
[[-12, 0, 0], [16, 1/3, 0], [49, -35/12, -1], [49, -35/12, 1]]
[[-12, 0, 0], [17, 5/12, 0], [41, -143/60, -2/5], [44, -32/15, 8/5]]
[[-12, 0, 0], [21, 3/4, 0], [28, -4/3, 0], [37, -7/12, 2]]
[[-12, 0, 0], [21, 3/4, 0], [29, -5/4, -2/3], [32, -1, 4/3]]
[[-12, 0, 0], [25, 13/12, 0], [25, -119/156, -10/13], [28, -20/39, 16/13]]
[[-13, 0, 0], [14, 1/13, 0], [182, -13, 0], [183, -168/13, 2]]
[[-13, 0, 0], [15, 2/13, 0], [98, -84/13, -1], [98, -84/13, 1]]
[[-13, 0, 0], [18, 5/13, 0], [47, -168/65, -2/5], [50, -153/65, 8/5]]
[[-13, 0, 0], [23, 10/13, 0], [30, -84/65, -1/5], [38, -44/65, 9/5]]
[[-14, 0, 0], [15, 1/14, 0], [210, -14, 0], [211, -195/14, 2]]
[[-14, 0, 0], [18, 2/7, 0], [63, -7/2, 0], [67, -45/14, 2]]
[[-14, 0, 0], [19, 5/14, 0], [54, -96/35, -4/5], [55, -187/70, 6/5]]
[[-14, 0, 0], [22, 4/7, 0], [39, -12/7, -1/2], [43, -10/7, 3/2]]
[[-14, 0, 0], [27, 13/14, 0], [31, -171/182, -10/13], [34, -66/91, 16/13]]
[[-15, 0, 0], [16, 1/15, 0], [240, -15, 0], [241, -224/15, 2]]
[[-15, 0, 0], [17, 2/15, 0], [128, -112/15, -1], [128, -112/15, 1]]
[[-15, 0, 0], [24, 3/5, 0], [40, -5/3, 0], [49, -16/15, 2]]
[[-15, 0, 0], [24, 3/5, 0], [41, -8/5, -2/3], [44, -7/5, 4/3]]
[[-15, 0, 0], [28, 13/15, 0], [33, -72/65, -6/13], [40, -25/39, 20/13]]
[[-15, 0, 0], [32, 17/15, 0], [32, -161/255, -16/17], [33, -48/85, 18/17]]

• Your example illustration seems to have only included the "inside" apollonian circles after the first operation. – Sparr Oct 8 '14 at 20:48
• @Sparr I'm not sure what you mean. After the first operation, one of the two Apollonian circles already exists (the original parent circle that you didn't pick for the current iteration) and you're only looking for the other solution. – Martin Ender Oct 8 '14 at 20:52
• Never mind, you're right, I was mis-reading. – Sparr Oct 8 '14 at 23:37

GolfScript (289 bytes vector / 237 bytes raster)

At 289 bytes and executing in a reasonable time:

'/'/n*','/']['*0,1/*~1.$[]*(~-400*:&;{1+1=*}/:D;{{1+2<~D@*\/}%}%'<svg><g fill="none" stroke="red">'puts.{[[~@:b[D&*\abs]{@&*[b]+}2*]{'.0/'*'"#{ }"'n/*~}%'<circle r=" " cx=" " cy=" " />'n/\]zip puts}:|/[{.([.;]+}3*]{(:?zip{)\~++2*\-}%:c.|0=D&*<{?);[c]+[{([.;]+.}3*;]+}*.}do'</g></svg>'  This takes input on stdin and generates an SVG file to stdout. Unfortunately it takes a bit too long for an online demo, but a tweaked version which aborts early can give you an idea. Given input [[-2, 0, 0], [3, 1/2, 0], [6, -2, 0], [7, -3/2, 2]] the output (converted to PNG with InkScape) is At 237 bytes and taking far too long (I extrapolate that it would take just over a week to produce similar output to the above, although in one-bit black and white): '/'/n*','/']['*0,1/*~1.$[]*(~-400*:&;{1+1=*}/:D;{{1+2<~D@*\/}%}%.[{.([.;]+}3*]{(:?[zip{)\~++2*\-}%:c]@+\0c=D&*<{?);[c]+[{([.;]+.}3*;]+}*.}do;:C;'P1 ''801 '2*.~:B*,{:P;C{:?[0=2/.D&*-.*\D&*+.*]{2,{P{B/}2$*B%400-?0=*\)?=&*-.*}/+<},,1=},!}/  Output is NetPBM format without newlines, so possibly doesn't strictly follow the spec, although the GIMP will still load it. If strict conformance is required, insert an n after the last !. The rasterisation is by testing each pixel against each circle, so the time taken is pretty much linear in the number of pixels times the number of circles. By downscaling everything by a factor of 10, '/'/n*','/']['*0,1/*~1.$[]*(~-40*:&;{1+1=*}/:D;{{1+2<~D@*\/}%}%.[{.([.;]+}3*]{(:?[zip{)\~++2*\-}%:c]@+\0c=D&*<{?);[c]+[{([.;]+.}3*;]+}*.}do;:C;'P1 ''81 '2*.~:B*,{:P;C{:?[0=2/.D&*-.*\D&*+.*]{2,{P{B/}2\$*B%40-?0=*\)?=&*-.*}/+<},,1=},!}/


will run in 10 minutes and produce

(converted to PNG with the GIMP). Given 36 hours it produced the 401x401

• I never would've thought you could do graphical output with Golfscript... – Beta Decay Oct 13 '14 at 17:49

JavaScript (418 410 bytes)

Implemented as a function:

function A(s){P='<svg><g fill=none stroke=red transform=translate(400,400)>';Q=[];s=eval(s);S=-400*s[0][0];function d(c){P+='<circle r='+Math.abs(p=S/c[0])+' cx='+p*c[1]+' cy='+p*c[2]+' />'}for(c=4;c--;d(s[0]),s.push(s.shift()))Q.push(s.slice());for(;s=Q.shift();d(c)){c=[];for(i=4;i--;)c[i]=2*(s[0][i]+s[1][i]+s[2][i])-s[3][i];for(i=6;c[0]<S&&i;)Q.push([s[i--%3],s[i--%3],c,s[i%3]])}document.body.innerHTML=P}


Online demo (note: doesn't work in browsers which fail to honour the SVG spec's requirements with respect to implicit sizing, so I offer a slightly longer version which works around that bug; browsers may also render the SVG less accurately than e.g. Inkscape, although Inkscape is a bit stricter on quoting attributes).

Note that 8 bytes could be saved by using document.write, but that seriously borks jsFiddle.

• You can probably save more by defining the function with ES6 and storing, e.g., S/c[0] in a variable and then also get rid of Math.abs with a ternary operator etc. – Ingo Bürk Oct 8 '14 at 18:19
• @IngoBürk, if I were going to go the ES6 route then I'd write it in CoffeeScript instead. – Peter Taylor Oct 8 '14 at 18:34
• use the host c99.nl. It allows document.write. – xem Oct 8 '14 at 19:51
• Good to see an answer to this:) – MickyT Oct 8 '14 at 20:04
• Updated with @IngoBürk's suggestion for a temporary variable. Eliminating Math.abs would actually cost a character. – Peter Taylor Oct 10 '14 at 15:56

Mathematica 289 characters

By solving the bilinear system as per http://arxiv.org/pdf/math/0101066v1.pdf Theorem 2.2 (highly inefficient).

Spaces not needed, still golfing it:

w = {k, x, y};
d = IdentityMatrix;
j = Join;
p_~f~h_ := If[#[[-1, 1]] < 6! h,
q = 2 d@4 - 1;
m = #~j~{w};
r = Complement[w /. NSolve[ And @@ j @@
If[r != {},
a~AppendTo~# & @@ r;
Function[x, x~j~{#}~f~h & /@ r]@#]] & /@ p~Subsets~{3};
Graphics[Circle @@@ ({{##2}, 1}/# & @@@ (f[a = #, -Tr@#]; a))] &


A reduced size animation with input {{-13, 0, 0}, {23, 10/13, 0}, {30, -84/65, -1/5}, {38, -44/65, 9/5}}

• How do you take input? – Martin Ender Oct 9 '14 at 14:26
• @MartinBüttner as a function argument, by adding @{{-1, 0, 0}, {2, 1, 0}, {2, -1, 0}, {3, 0, 2}} to the last line – Dr. belisarius Oct 9 '14 at 14:27
• @MartinBüttner If you're going to test it try first with 50/h instead of 400/h. You're going to get the result faster. also, you can monitor the progress by entering Dynamic@Length@a before executing the function – Dr. belisarius Oct 9 '14 at 14:29
• Instructions for testing this answer (with a reduced number of circles) without Mathematica installed: 1) Download this from pastebin and save it as *.CDF 2) Download and install the free CDF environment from Wolfram Research at (not a small file). Enjoy. Tell me if it works!-- Note: Calcs are slow, wait for the graphics to come up. – Dr. belisarius Oct 9 '14 at 22:47
• What does the "highly inefficient" comment refer to? Is it that (looking at the animation) you're apparently drawing most of the circles at least twice? I think the complex Descartes approach is inherently as efficient as it gets. – Peter Taylor Oct 10 '14 at 12:00

Maple (960 bytes)

I used Descartes Theorem to generate the Apollonian Gasket and then use Maple's plotting system to plot it. If I have time I want to further golf this and and change it to Python (Maple is definitely not the best for fractals). Here is a link to a free Maple player if you want to run my code.

X,Y,Z,S,N:=abs,evalf,member,sqrt,numelems;
f:=proc(J)
L:=map((x)->[x[1],(x[2]+x[3]*I)/x[1]+50*(1+I)/X(J[1][2])],J);
R:=Vector([L]);
T,r:=X(L[1][3]),L[1][4];
A(L[1][5],L[2][6],L[3][7],L[1][8],L[2][9],L[3][10],R,T,r);
A(L[1][11],L[2][12],L[4][13],L[1][14],L[2][15],L[4][16],R,T,r);
A(L[1][17],L[3][18],L[4][19],L[1][20],L[3][21],L[4][22],R,T,r);
A(L[2][23],L[3][24],L[4][25],L[2][26],L[3][27],L[4][28],R,T,r);
plots[display](seq(plottools[circle]([Re(R[i][29]),Im(R[i][30])],X(1/R[i][31])),i=1..N(R))):
end proc:
A:=proc(a,b,c,i,j,k,R,E,F)
K:=i+k+j+2*S(i*k+i*j+k*j);
if K>400*E then
return;
end if;
C:=(a*i+c*k+b*j+2*S(a*c*i*k+b*c*j*k+a*b*i*j))/K;
C2:=(a*i+c*k+b*j-2*S(a*c*i*k+b*c*j*k+a*b*i*j))/K;
if Y(X(C-F))<1/E and not Z([K,C],R) then
R(N(R)+1):=[K,C];
A(a,b,C,i,j,K,R,E,F);
A(a,c,C,i,k,K,R,E,F);
A(b,c,C,j,k,K,R,E,F);
end if:
if Y(X(C2-F))<1/E and not Z([K,C2],R) then
R(N(R)+1):=[K,C2];
A(a,b,C2,i,j,K,R,E,F);
A(a,c,C2,i,k,K,R,E,F);
A(b,c,C2,j,k,K,R,E,F);
end if:
end proc:


f([[-1, 0, 0], [2, 1, 0], [2, -1, 0], [3, 0, 2]]);

f([[-9, 0, 0], [14, 5/9, 0], [26, -77/45, -4/5], [27, -8/5, 6/5]]);
`