Fractals/Apollonian fractals

From Wikibooks, open books for an open world
Jump to: navigation, search

Theory[edit]

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

How to compute and draw Apollonian gasket ?[edit]

Notation[edit]

  • 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[edit]

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 aproximation of limit set which is made of infinite number of stages (cicles).

Mandelbrot algorithm[edit]

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

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 \pm2 \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_1c_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[edit]

Functions for computing curvature and center[edit]
 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[edit]

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

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

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 equall 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_0out - r_i)*i

c_{ib} = (r_{0out} + r_i) + (2*r_0out - r_i - heght)*i

c_{ic} = (r_{0out} - r_i) + (2*r_0out - r_i - heght)*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 betwen circles[edit]

 {\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[edit]

  • 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 cofiguration ( 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[edit]

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

Java[edit]

CirclePack program by Ken Stephenson

Maxima CAS[edit]

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

Common Lisp[edit]

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[edit]
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()
]]></script></svg>

References[edit]

  1. Apollonian gasket
  2. mathforum : Problem of Apollonius
  3. Apollonius' Tangency Problem at MathPages
  4. Mandelbrot's Algorithm at Fractal Geometry from Yale University by Michael Frame, Benoit Mandelbrot, and Nial Neger
  5. Drawing 2D Apollonian gasket with n identical circles in Matlab by Guillaume JACQUENOT
  6. Apollonian Gasket in Mathematica and Povray by Paul Nylander
  7. Drawing the Apollonian Gasket with Common Lisp and Vecto by Luis Diego Fallas
  8. Leibnitz packing by Takaya Iwamoto with program in AutoLisp
  9. Apollonian Gasket - by Paul Bourke in BAsic and C
  10. How to pack n circles inside unit circle by Erich Friedman
  11. SVG Math Animation Example: Apollonian Gasket (-15,32,32,33) by KiHyuck Hong. See the src code.
  12. SVG Math Animation Example: Apollonian Gasket (-15,32,32,33) by KiHyuck Hong
  13. File Apolonian gasket containes about 8140 circles, but it looks like made of about 20 steps. So it should have about 10 460 353 205 circles