Theory

wikipedia^[1]
mathforum : Problem of Apollonius ^[2]

How to compute and draw Apollonian gasket ?

Notation

Apollonian gasket = curvilinear Sierpinski gasket = Leibniz packing = Apollonian packing
stage = level = step
incircle = inscribed circle = circle inside 3 circles = inner circle
gap = curvilinear triangle = ideal triangle ( because the sides are tangent at each vertix and the angle between them is zero) = region between 3 tangent circles

Algorithms

Ancillary figure for van Roomen's solution to the problem of Apollonius.

Apollonian gasket can be made using these algorithms:^[3]

Mandelbrot algorithm ^[4]( using circle inversion)^[5]^[6]
Soddy's algorithm ( using circle Descartes theorem )^[7]^[8]
IFS algorithm by Alexei Kravchenko and Dmitriy Mekhontsev ^[9]

All algorithms give state after n stages. It is an approximation of limit set which is made of infinite number of stages (cicles).

Mandelbrot algorithm

Gasket made from 4 equal inner circles made with circle inversion

Animation of the solution to Apollonius' problem using the method of inversion.

(to do )

Soddy's algorithm

Integer Apollonian gasket defined by circle curvatures of (32,32,33). This image is probably made of 19 stages

It will be explained by example gasket with initial configuration : 3 inner circles with integer curvatures

(k_{ia},k_{ib},k_{ic})=(32,32,33)

This algorithm uses:

previous circles defined by centers $c_{n}$ and curvatures $k_{n}$
circle Descartes theorem to compute curvature of next circle $k_{4}$

k_{4}=k_{1}+k_{2}+k_{3}\pm 2{\sqrt {k_{1}k_{2}+k_{2}k_{3}+k_{3}k_{1}}}

complex circle Descartes theorem to compute center of next circle $c_{4}$

c_{4}={\frac {c_{1}k_{1}+c_{2}k_{2}+c_{3}k_{3}\pm 2{\sqrt {k_{1}k_{2}c_{1}c_{2}+k_{2}k_{3}c_{2}c_{3}+k_{1}k_{3}c_{1}c_{3}}}}{k_{4}}}

When center of outer circle is at origin one has to chose solutions in case of centers. It is easier when the all gasket is in the first quadrant of Cartesian plane. Then all centers have positive parts ( real and imaginary ).

Functions

Functions for computing curvature and center

f_{p}(a_{1},a_{2},a_{3}){\overset {\underset {\mathrm {def} }{}}{=}}a_{1}+a_{2}+a_{3}+2\cdot {\sqrt {a_{1}a_{2}+a_{2}a_{3}+a_{3}a_{1}}}

f_{\color {Red}m}(a_{1},a_{2},a_{3}){\overset {\underset {\mathrm {def} }{}}{=}}a_{1}+a_{2}+a_{3}{\color {Red}-}2\cdot {\sqrt {a_{1}a_{2}+a_{2}a_{3}+a_{3}a_{1}}}

Note that variable $a_{i}$ can be curvature $k_{i}$ or $c_{i}*k_{i}$ so $f_{p}$ can be used for computing both curvatures and centers

Functions for computing ck list

It is easier to use a list $ck$ which defines circle. It consist of two element center $c$ and curvature $k$ :

$ck=[c,k]$

$c=ck[1]$

$k=ck[2]$

Now one can define new functions ( Maxima CAS code) for computing 4-th circle mutually tangent to 3 given circles:

f_pp(ck1,ck2,ck3):=block
([c4,k4,ck4],
k4:f_p(ck1[2],ck2[2],ck3[2]),
c4:f_p(ck1[1]*ck1[2],ck2[1]*ck2[2],ck3[1]*ck3[2])/k4,
ck4:[c4,k4])$

f_pm(ck1,ck2,ck3):=block
([c4,k4,ck4],
 k4:f_p(ck1[2],ck2[2],ck3[2]),
 c4:f_m(ck1[1]*ck1[2],ck2[1]*ck2[2],ck3[1]*ck3[2])/k4,
 ck4:[c4,k4]);

f_mm(ck1,ck2,ck3):=block
([c4,k4,ck4],
 k4:f_m(ck1[2],ck2[2],ck3[2]),
 c4:f_m(ck1[1]*ck1[2],ck2[1]*ck2[2],ck3[1]*ck3[2])/k4,
 ck4:[c4,k4]);

 f_mp_ck4(ck1,ck2,ck3):=block
([c4,k4],
k4:solve_m_eq(ck1[2],ck2[2],ck3[2]),
c4:solve_p_eq(ck1[1]*ck1[2],ck2[1]*ck2[2],ck3[1]*ck3[2])/k4,
return([c4,k4])
);

Function $f_{pp}$ uses $f_{p}$ function for computing both center and curvature. It is used for computing almost all circles without outer circle and hard circles

Function $f_{pm}$ uses $f_{p}$ function for computing curvature and $f_{m}$ function for computing center. It is used for computing hard circles.

Function $f_{mm}$ uses $f_{m}$ function for computing both center and curvature. It is used for computing outer circle.

Functions for filling gaps

fill_easy(ckla,cklb,cklc,max_stage):=block
([ckm],
 if max_step>0 then (
  /* circle in the middle */
  ckm:give_pp_ck4(ckla,cklb,cklc),
  /* 3 circles around */
  fill_easy(ckla,cklb,ckm,max_stage-1),
  fill_easy(cklb,cklc,ckm,max_stage-1),
  fill_easy(ckla,cklc,ckm,max_stage-1)))$

fill_hard(ck_ic,ck_ia,ck_0out,max_stage):=block
([ck_1c,ck_2cc,ck_3ccc,ck_4ccca,ck_5cccac], /* hard circles */
 if 	max_stage>0 then (
  /* step 1 = circle in the middle */
  ck_1c:f_pm(ck_ia,ck_ic,ck_0out),      /* 1c */
  /* step 2 = 3 circles around */
  fill_easy(ck_ia,ck_0out,ck_1c,max_stage-1), /* 2ca */
  fill_easy(ck_ia,ck_ic,ck_1c,max_stage-1), /* 2cb */
  /* hard subgap 2cc */
  if max_stage>1 then ck_2cc:f_pm(ck_ic,ck_1c,ck_0out), /* 2cc */
  /* step 3 for subgap 2cc */
  fill_easy(ck_0out,ck_1c,ck_2cc,max_stage-2), /* 3cca */
  fill_easy(ck_ic,ck_1c,ck_2cc,max_stage-2),  /* 3ccb */
  if max_stage>2 then ck_3ccc:f_pm(ck_ic,ck_0out,ck_2cc),	/* 3ccc */
  /* step 4 for subgap 3ccc */
  if max_stage>3 then ck_4ccca:f_pm(ck_0out,ck_2cc,ck_3ccc),
  fill_easy(ck_ic,ck_2cc,ck_3ccc,max_stage-3),  /* 4cccb */
  fill_easy(ck_ic,ck_0out,ck_3ccc,max_stage-3),  /* 4cccc */
  /* step 5 for subgap 4ccca */
  fill_easy(ck_0out,ck_2cc,ck_4ccca,max_stage-4),  /*  5cccaa */
  fill_easy(ck_2cc,ck_3ccc,ck_4ccca,max_stage-4),  /* 5cccab */
  if max_stage>4 then ck_5cccac:f_pm(ck_0out,ck_3ccc,ck_4ccca)))$

Construction

Apollonian gasket with marked hard circles

initial stage : 3 mutually tangent circles = $(k_{ia},k_{ib},k_{ic})$
zero stage : 2 circles tangent to previous : outer circle $ck_{0out}$ and inner circle $ck_{0in}$ , where :

{\color {Red}ck}_{0out}=f_{\color {Red}mm}(ck_{ia},ck_{ib},ck_{ic})

ck_{0in}=f_{pp}(ck_{ia},ck_{ib},ck_{ic})

first stage : there are 6 gaps ( = curvilinear triangles) between previous circles. Three gaps are peripheral ( near outer circle) and three are medial ( near inner circle). Put one circle in each gap :

ck_{1a}=f_{pp}(ck_{ia},ck_{ib},ck_{0out})

ck_{1b}=f_{pp}(ck_{ib},ck_{ic},ck_{0out})

{\color {Red}c}k_{1c}=f_{p{\color {Red}m}}(ck_{ic},ck_{ia},ck_{0out})

ck_{1d}=f_{pp}(ck_{ia},ck_{ib},ck_{0in})

ck_{1e}=f_{pp}(ck_{ib},ck_{ic},ck_{0in})

ck_{1f}=f_{pp}(ck_{ic},ck_{ia},ck_{0in})

Each new circle placed in the gap makes three new gaps (trifurcation).

second stage : there are 18 (6*3) new gaps. Put one circle in each gap using $f_{pp}$ function except one hard circle :

{\color {Red}c}k_{2cc}:f_{p{\color {Red}m}}(ck_{ic},ck_{1c},ck_{0out})

next stages : Put one circle in each gap at each stage. At each stage there is only one hard circle and it is in $1c$ gap.

Algorithm for programmers

Because all circles are inscribed in outer circle so it is easier to start from it.

choose size of quadratic image $side$
compute radius of outer circle :

r_{0out}={\frac {side}{2}}

place outer circle in first quadrant of Cartesian plane :

c_{0out}=r_{0out}+r_{0out}*i

compute radius of 3 equal initial circles ^[10]

a=1+{\frac {2}{\sqrt {3}}}

r_{i}={\frac {r_{0out}}{a}}

find centers of initial circles. For example place one circle $ia$ in upper row and 2 circles ( $ib$ and $ic$ ) in lower row.

Then centers of these circles are vertices of equilateral triangle, with side of the triangle $=2*r_{i}$ so :

$height=r_{i}*{\sqrt {3}}$

$c_{ia}=r_{0out}+(2*r_{0}out-r_{i})*i$

$c_{ib}=(r_{0out}+r_{i})+(2*r_{0}out-r_{i}-height)*i$

$c_{ic}=(r_{0out}-r_{i})+(2*r_{0}out-r_{i}-height)*i$

compute inner circle $0in$

ck_{0in}=f_{pp}(ck_{ia},ck_{ib},ck_{ic})

"Now we have a configuration to start" ^[11] There is 6 gaps to fill. 5 gaps are easy and one gap $1c$ is hard to fill.

Relation between circles

{\color {Red}ck}_{0out}=f_{\color {Red}mm}(ck_{ia},ck_{ib},ck_{ic})

ck_{0in}=f_{pp}(ck_{ia},ck_{ib},ck_{ic})

{\color {Red}ck}_{ia}=f_{\color {Red}mm}(ck_{ib},ck_{ib},ck_{0out})

c{\color {Red}k}_{ib}=f_{{\color {Red}m}p}(ck_{ia},ck_{ic},ck_{0out})

{\color {Red}ck}_{ic}=f_{\color {Red}mm}(ck_{ia},ck_{ib},ck_{0out})

Number of circles

number of new circles at stage n = $2*3^{n}$
total number of circles after n stages = $2+3^{n+1}$

where $n>=0$

For example :

stage -1 = initial configuration ( 3 inner circles )
stage 0 gives 2 new circles ( one outer and one most inner ). All circles = 5
stage 1 gives 6 new circles ( all circles = 11 )
stage 2 gives 18 new circles ( all circles = 29 )
stage 3 gives 54 new circles ( all circles = 83)
stage 4 gives 162 new circles ( all circles = 245)
stage 5 gives 486 new circles ( all circles = 731)
stage 6 gives 1 458 new circles ( all circles = 2 189 )
stage 7 gives 4 374 new circles ( all circles = 6 563 )
stage 8 gives 13 122 new circles ( all circles = 19 685 )
stage 9 gives 39 366 new circles ( all circles = 59 051 )
stage 10 gives 118 098 new circles ( all circles = 177 149 )
stage 11 gives 354 294 new circles ( all circles = 531 443 )
stage 12 gives 1 062 882 new circles ( all circles = 1 594 325 )
stage 13 gives 3 188 646 new circles ( all circles = 4 782 971 )
stage 14 gives 9 565 938 new circles ( all circles = 14 348 909 )
stage 15 gives 28 697 814 new circles ( all circles = 43 046 723 )
stage 16 gives 86 093 442 new circles ( all circles = 129 140 165 )
stage 17 gives 258 280 326 new circles ( all circles = 387 420 491 )
stage 18 gives 774 840 978 new circles ( all circles = 1 162 261 469 )
stage 19 gives 2 324 522 934 new circles ( all circles = 3 486 784 403 )
stage 20 gives 6 973 568 802 new circles ( all circles = 10 460 353 205 )

It can be computed using this Maxima CAS code :

 NumberOfStageCircles(stage):= 
  if stage =-1 then 3
  elseif stage=0 then 2
  else 2*3^stage;

 NumberOfAllCircles(stage_max):=sum(NumberOfStageCircles(i),i,-1,stage_max);

for i:-1 step 1 thru 20 do print("stage ",i," gives ",NumberOfStageCircles(i)," new circles ( all circles =  ",NumberOfAllCircles(i), " )");

Number of circles grows rapidly and makes problems caused by big file size. One can :

use only stages up to 5-7. It seems to be sufficient to give good image in a short time.^[12]
do not draw circles which are small ( have a small radius < pixel size) ^[13]
convert svg file to png ( for example using Image Magic : "convert a.svg a.png" ). Then for stage 12 svg image has 128 MB and png file only 57 KB.
be happy that you have so detailed picture and wait a long time to load file or see how system hangs up (:-))

Curvatures

Note that only one curvature has been computed with $f_{\color {Red}m}$ function.

It is $k_{0out}$ . It is also only one negative curvature.

All the rest have been comuted with $f_{p}$ function.

stage $n\,$	curvatures
-1	$(32,32,33)\,$
0	$({\color {Red}-15},209)\,$
1	$(68,65,68,516,513,516)\,$
2	$(953,1481,1484,945,1476,1476,953,1484,1481,281,137,140,273,132,132,68,137,140)\,$

Example programs

Haskell and EDSL Diagrams

Image and code by Brent Yorgey^[14]
diagrams code ^[15]

Java

CirclePack program by Ken Stephenson

Maxima CAS

Program is explained above, and src is in description of Apollonian_rc_6.svg

Common Lisp

To run put apollonian.lisp file in your home directory and in Emacs/slime. To get another step change *max-step*

Load file to run code ( batch mode ) :

(load  "apollonian.lisp")

; first run
;(require :asdf)
;(require :asdf-install)
;(asdf-install:install :zpng)
;(asdf-install:install :Vecto)
; you must press 2 and 0 when the program asks

; next times load packages from disk
(asdf:operate 'asdf:load-op 'vecto)

(in-package vecto)

;-----------------------------------------------------------------------
;---------------------------- definitions of functions ---------------
;-----------------------------------------------------------------------

; kc is a list describing circle (list k c)  where k is curvature and c is a center ( complex number)
; k can be positive ( inner circle) , negative ( outer circle) or zero ( =  line )

(defun draw-vecto-circle (kc)
  (let* ((c (second kc))
	 (r (abs (/ 1 (first kc)))))	; r = abs(1/k)
    (centered-circle-path (realpart c) (imagpart c) r))) ; vecto

; Desfirsttes' circle theorem 
(defun solve-equation-m (k1 k2 k3)
  (- (+ k1 k2 k3 ) (* 2 (sqrt (+ (* k1 k2) (* k2 k3) (* k1 k3))))))

(defun solve-equation-p (k1 k2 k3)
  (+ (+ k1 k2 k3 ) (* 2 (sqrt (+ (* k1 k2) (* k2 k3) (* k1 k3))))))

; compute, draws and gives back 4-th circle usind solve-equation-p and solve-equation-p
(defun draw-forth-circle-pp (kc1 kc2 kc3)
  (let* ((k1 (first kc1))
	 (k2 (first kc2))
	 (k3 (first kc3))
	 (p1 (* k1 (second kc1))) ; pn = kn * cn 
	 (p2 (* k2 (second kc2)))
	 (p3 (* k3 (second kc3)))
	 (k4 (solve-equation-p k1 k2 k3)) ; find k
	 (c4 (/ (solve-equation-p p1 p2 p3) k4)) ; find c
	 (kc4 (list k4 c4)))
    (draw-vecto-circle kc4) ; draw 
    kc4)) ; return 4-th circle

; compute, draws and gives back 4-th circle usind solve-equation-p and solve-equation-m
(defun draw-forth-circle-pm (kc1 kc2 kc3)
  (let* ((k1 (first kc1))
	 (k2 (first kc2)) 
	 (k3 (first kc3))
	 (p1 (* k1 (second kc1)))
	 (p2 (* k2 (second kc2)))
	 (p3 (* k3 (second kc3)))
	 (k4   (solve-equation-p k1 k2 k3))
	 (c4 (/ (solve-equation-m p1 p2 p3) k4))
	 (kc4 (list k4 c4)))
    (draw-vecto-circle kc4)
    kc4))

; easy = pp
(defun fill-easy-gap (kc1 kc2 kc3 steps)
  (when (> steps 0)
      (let* ((kc4 (draw-forth-circle-pp kc1 kc2 kc3))) ; 4-th circle
	(fill-easy-gap kc1 kc2 kc4 (- steps 1)) ; 3 subgaps
	(fill-easy-gap kc2 kc3 kc4 (- steps 1))  
	(fill-easy-gap kc3 kc1 kc4 (- steps 1)))))

; only for gap 1c = hard gap
(defun fill-hard-gap (kc-ic kc-ia kc-0out steps)
  (let* (kc-1c kc-2cc  kc-3ccc  kc-4ccca kc-5cccac) ; hard circles , also kc-6cccacc
    ;------------ gap 1c --------------------------------------------   
    (when (> steps 0) (setq kc-1c (draw-forth-circle-pm kc-ic kc-ia kc-0out)) 						; kc-1c hard
	  (fill-easy-gap kc-ia kc-0out kc-1c (- steps 1))								; kc-2ca
	  (fill-easy-gap kc-ia kc-ic kc-1c (- steps 1))									; kc-2cb
	  ; ----------- hard subgap 2cc ------------------------------------
	  (when (> steps 1)  (setq kc-2cc (draw-forth-circle-pm kc-ic kc-0out kc-1c))					; kc-2cc hard
		(fill-easy-gap kc-0out kc-1c kc-2cc (- steps 2))							; kc-3cca
		(fill-easy-gap kc-ic kc-1c kc-2cc (- steps 2))	        						; kc-3ccb
		; ----------- hard subgap 3ccc ------------------------------------
		(when (> steps 2)   (setq kc-3ccc (draw-forth-circle-pm kc-ic kc-0out kc-2cc))				; kc-3ccc hard
		      (fill-easy-gap kc-ic kc-2cc kc-3ccc (- steps 3))							; kc-4cccb
		      (fill-easy-gap kc-ic kc-0out kc-3ccc (- steps 3))							; kc-4cccc
		      ; ----------- hard subgap 4ccca ------------------------------------
		      (when (> steps 3)  (setq kc-4ccca (draw-forth-circle-pm kc-0out kc-2cc kc-3ccc))			; kc-4ccca hard
			    (fill-easy-gap kc-0out kc-2cc kc-4ccca (- steps 4))						; kc-5cccaa
			    (fill-easy-gap kc-2cc kc-3ccc kc-4ccca (- steps 4))						; kc-5cccab
			    (when (> steps 4) (setq kc-5cccac (draw-forth-circle-pm kc-0out kc-3ccc kc-4ccca)))))))))   ; kc-5cccac hard

; example of use : (draw-apollonian-gasket "a.png" 800 2)
; classical Desfirsttes configuration : 3 mutually tangent circles and 4-th inside ( and fith outside)
(defun draw-apollonian-gasket-n3 (file side step)
  (with-canvas (:width side :height side) ; vecto
	       (set-rgb-stroke 0 0 0) ; vecto
	       (set-line-width 1) ; vecto

	       
	      (let* (	(r0out (/ side 2))
			(k0out (- (/ 1 r0out)))
			(c0out (complex r0out r0out))
			(kc0out (list k0out c0out)) ; list defining circle 
					; a1:float(1 + 2 / sqrt(3)),  http://www2.stetson.edu/~efriedma/cirincir/ 
			(a (+ 1 (/ 2 (sqrt 3))))
			(ri   (/ r0out a)) ; three equal inner circles = initial circles
			(ki (/ 1 ri))	
			(cia (complex r0out (- side ri))) ; one circle in upper row
			(kcia (list ki cia))
			(h (* ri (sqrt 3))) ; height of equilateral triangle with side = 2*ri
					; 2 circles in lower row
			(cib (complex  (+ r0out ri) (- side (+ h ri))))
			(cic (complex  (- r0out ri) (- side (+ h ri))))
			(kcib (list ki cib))
			(kcic (list ki cic))
			kc0in ) 
			; drawing code
			; step -1 = three equal inner circles
		 (draw-vecto-circle kcia ) 		; ia
		 (draw-vecto-circle kcib ) 		; ib
		 (draw-vecto-circle kcic ) 		; ic
			; step 0 = outer and most inner circle
		 (draw-vecto-circle kc0out )		; 0out
		 (setq kc0in ( draw-forth-circle-pp kcia kcib kcic )) ; 0in
			; this is starting configuration. 
			; One can go to the next steps now
			; Fill 6 gaps :   
			; 3 peripheral gaps
		 (fill-easy-gap kcia kcib kc0out step) ; 1a
		 (fill-easy-gap kcib kcic kc0out step) ; 1b	
		 (fill-hard-gap kcic kcia kc0out step) ; 1c
			; 3 medial gaps
		 (fill-easy-gap kcia kcib kc0in step)  ; 1d
		 (fill-easy-gap kcib kcic kc0in step)  ; 1e
		 (fill-easy-gap kcic kcia kc0in step)  ; 1f
			; rest of drawing code
		 (stroke)) ; before save ( vecto procedure)
	       (print (save-png file)))) ; info text and vecto procedure

;---------------compile ---------------------------------
(compile 'draw-vecto-circle )
(compile 'solve-equation-m)
(compile 'solve-equation-p)
(compile 'draw-forth-circle-pp)
(compile 'draw-forth-circle-pm)
(compile 'fill-easy-gap)
(compile 'fill-hard-gap)
(compile 'draw-apollonian-gasket-n3)

;----------global var ----------------------

(defparameter *max-step* 0 " maximal step of drawing apollonian gasket from. It is an integer >= 0 ")

(defparameter *file-name*
  (make-pathname 
   :name (concatenate 'string "apollonian-n3-s" (write-to-string *max-step*))
   :type "png")
  "name (or pathname) of png file ")

;====================================== run ==================================
(draw-apollonian-gasket-n3 *file-name* 800 *max-step*)

JavaScript

Random Apollonian circle fractal

This is an incremental algorithm that maintains a queue of triples of touching circles. The main loop starts at the front of the queue, removing a triple and filling it with the required circle, adding the resulting circle to a list of circles to draw as SVG and adding the three triples that result from adding the circle to the back of the computation queue. As a bonus feature, the circles are coloured so that no two neighbours have the same colour.

Usage: concatenate the two SVG boilerplate parts and the middle Javascript part into one SVG file and load into a web browser (or another SVG viewer that supports Javascript or ECMAScript).

<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
<svg
 xmlns="http://www.w3.org/2000/svg"
 xmlns:xlink="http://www.w3.org/1999/xlink"
 width="100%" height="100%" viewBox="0 0 768 768"
 onload="init(evt)"
><title>Apollonian Fractal</title
><desc>Randomized Apollonian fractal circle packing</desc
><script type="text/ecmascript"><![CDATA[

// SVG properties
var ns = "http://www.w3.org/2000/svg";
var xns = "http://www.w3.org/1999/xlink";
var SVGDocument = null;
var size = 768; // should match 'viewBox'
var minsize = 0.001;

// convert internal coordinates and lengths to SVG
function coord(x) { return x*size/2 + size/2; }
function dist(d) { return d*size/2; }

// complex numbers
function C(r,i)    { this.r = r; this.i = i; }
function Cconj(c)  { return new C(c.r, -c.i); }
function Cadd(c,d) { return new C(c.r+d.r,c.i+d.i); }
function Csub(c,d) { return new C(c.r-d.r,c.i-d.i); }
function Cmul(c,d) { return new C(c.r*d.r-c.i*d.i,c.r*d.i+c.i*d.r); }
function Csqrt(c)  { return new C(Math.sqrt(Math.sqrt(c.r*c.r+c.i*c.i)) * Math.cos(Math.atan2(c.i,c.r)/2),
                                  Math.sqrt(Math.sqrt(c.r*c.r+c.i*c.i)) * Math.sin(Math.atan2(c.i,c.r)/2)); }

// sort points into order
function clockwise(a, b, c) {
  var x1 = a.x-b.x; var y1 = a.y-b.y;
  var x2 = c.x-b.x; var y2 = c.y-b.y;
  if ((x1*y2-y1*x2) >= 0) return [a,b,c]; else return [a,c,b];
}

// compute the 4th circle touching 3 circles, each of which touch the other two
function Kiss(a, b, c, initial) {
  var k1 = 1 / a.r; var z1 = new C(a.x, a.y); var kz1 = Cmul(new C(k1,0),z1);
  var k2 = 1 / b.r; var z2 = new C(b.x, b.y); var kz2 = Cmul(new C(k2,0),z2);
  var k3 = 1 / c.r; var z3 = new C(c.x, c.y); var kz3 = Cmul(new C(k3,0),z3);
  var k4p = k1 + k2 + k3 + 2*Math.sqrt(k1*k2 + k2*k3 + k3*k1);
  var k4m = k1 + k2 + k3 - 2*Math.sqrt(k1*k2 + k2*k3 + k3*k1);
  var kz4p = Cadd(Cadd(Cadd(kz1,kz2),kz3),Cmul(new C(2,0),
    Csqrt(Cadd(Cadd(Cmul(kz1,kz2),Cmul(kz2,kz3)),Cmul(kz3,kz1)))));
  var kz4m = Csub(Cadd(Cadd(kz1,kz2),kz3),Cmul(new C(2,0),
    Csqrt(Cadd(Cadd(Cmul(kz1,kz2),Cmul(kz2,kz3)),Cmul(kz3,kz1)))));
  var k4;
  var kz4;
  var k4b;
  var kz4b;
  var cs = new Array();
  if (k4p > k4m) { k4 = k4p; kz4 = kz4p; k4b = k4m; kz4b = kz4m; }
  else           { k4 = k4m; kz4 = kz4m; k4b = k4p; kz4b = kz4p; }
  var cc = new Circle(kz4.r/k4,kz4.i/k4,Math.abs(1/k4), 6 - a.colour - b.colour - c.colour);
  var dx = a.x - cc.x
  var dy = a.y - cc.y
  var dr = a.r + cc.r
  if (Math.abs(dx * dx + dy *dy - dr * dr) > 0.0001) {
    cc = new Circle(kz4b.r/k4,kz4b.i/k4,Math.abs(1/k4), 6 - a.colour - b.colour - c.colour);
  }
  cs[cs.length] = cc;
  if (initial) {
    cc = new Circle(kz4b.r/k4b,kz4b.i/k4b,Math.abs(1/k4b), 6 - a.colour - b.colour - c.colour);
    cs[cs.length] = cc;
  }
  return cs;
}

// called once on load
function init(evt) {
  // get document
  SVGDocument = evt.target.ownerDocument;

  // initial bounding circle
  var b = new Circle(0, 0, -1, 0);

  // insert two randomly positioned touching circles
  var tr = 1-Math.random()/2;
  var pa = Math.PI/2 - Math.asin(Math.random()*(1-tr)/tr);
  var px = tr * Math.cos(pa);
  var py = tr * Math.sin(pa);
  var pr = 1 - tr;
  var qr = (1 - pr) * (1 - Math.cos(pa + Math.PI/2))
         / (1 + pr - (1 - pr) * Math.cos(pa + Math.PI/2));
  var qx = 0;
  var qy = qr - 1;
  var p = new Circle(px, py, pr, 1);
  var q = new Circle(qx, qy, qr, 2);

  // a queue to contain triples of touching circles
  var queue = new Array();
  var cs = Kiss(b,p,q, true);
  queue[queue.length] = [b,p,cs[0]];
  queue[queue.length] = [b,cs[0],q];
  queue[queue.length] = [cs[0],p,q];
  queue[queue.length] = [b,p,cs[1]];
  queue[queue.length] = [b,cs[1],q];
  queue[queue.length] = [cs[1],p,q];

  // a queue to contain circles to draw
  var draw = new Array();
  draw[draw.length] = b;
  draw[draw.length] = p;
  draw[draw.length] = q;
  draw[draw.length] = cs[0];
  draw[draw.length] = cs[1];

  // add 10000 more circles to the draw queue
  // adding new triples to the compute queue
  var c;
  for (c = 0; c < Math.min(queue.length, 10000); c = c + 1) {
    var c0 = queue[c][0];
    var c1 = queue[c][1];
    var c2 = queue[c][2];
    var ncs = Kiss(c0,c1,c2);
    var nc = ncs[0];
    if (nc.r > minsize) {
      queue[queue.length] = [nc,c1,c2];
      queue[queue.length] = [c0,nc,c2];
      queue[queue.length] = [c0,c1,nc];
      draw[draw.length] = nc;
    }
  }

  // add all generated circles to an SVG <g> element
  var g = SVGDocument.createElementNS(ns, "g");
  g.setAttributeNS(null, "stroke", "black");
  g.setAttributeNS(null, "stroke-width", "1");
  var i;
  for (i = 0; i < draw.length; i = i + 1) {
    g.appendChild(draw[i].draw());
  }

  // link the <g>s into the DOM so they are displayed
  var svg = SVGDocument.documentElement;
  svg.appendChild(g);

} // init()

// circle class
function Circle(x, y, r, colour) {
  // properties of a circle
  this.x = x;
  this.y = y;
  this.r = r;
  this.colour = colour;
  // convert to svg
  this.draw = function() {
    var colours = ["white", "red", "yellow", "cyan"];
    var c = SVGDocument.createElementNS(ns, "circle");
    c.setAttributeNS(null, "fill", colours[this.colour]);
    c.setAttributeNS(null, "cx", coord(this.x));
    c.setAttributeNS(null, "cy", coord(this.y));
    c.setAttributeNS(null, "r",  dist(this.r));
    return c;
  };
} // Circle()