理論

維基百科^[1]
數學論壇：阿波羅尼奧斯問題 ^[2]

如何計算和繪製阿波羅尼奧斯墊片？

符號

阿波羅尼奧斯墊片 = 曲線型謝爾賓斯基墊片 = 萊布尼茨填充 = 阿波羅尼奧斯填充
階段 = 層級 = 步
內切圓 = 內接圓 = 3 個圓內的圓 = 內圓
間隙 = 曲線型三角形 = 理想三角形（因為邊在每個頂點處相切，它們之間的角度為零）= 3 個相切圓之間的區域

演算法

阿波羅尼奧斯墊片可以使用以下算法制作：^[3]

曼德爾布羅特演算法 ^[4]（使用圓反演）^[5]^[6]
索迪演算法（使用圓的笛卡爾定理）^[7]^[8]
阿列克謝·克拉夫琴科和德米特里·梅洪特塞夫的 IFS 演算法 ^[9]

所有演算法在 n 個階段後給出狀態。它是一個無限多個階段（圓圈）的極限集的近似值。

曼德爾布羅特演算法

（待辦事項）

索迪演算法

它將透過具有初始配置的示例墊片進行解釋：3 個具有整數曲率的內圓

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

該演算法使用

由中心 $c_{n}$ 和曲率 $k_{n}$ 定義的先前圓
圓的笛卡爾定理來計算下一個圓的曲率 $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}}}

復圓的笛卡爾定理來計算下一個圓的中心 $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}}}

當外圓圓心位於原點時，需要在圓心情況下選擇解。當整個墊片位於笛卡爾平面第一象限時，更容易。然後所有圓心都有正部分（實部和虛部）。

函式

用於計算曲率和圓心的函式

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

注意，變數 $a_{i}$ 可以是曲率 $k_{i}$ 或 $c_{i}*k_{i}$ 因此 $f_{p}$ 可用於計算曲率和圓心

用於計算 ck 列表的函式

使用列表 $ck$ 更容易定義圓。它包含兩個元素：圓心 $c$ 和曲率 $k$

$ck=[c,k]$

$c=ck[1]$

$k=ck[2]$

現在可以定義新的函式（Maxima CAS 程式碼）來計算與給定的 3 個圓相切的第 4 個圓。

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

函式 $f_{pp}$ 使用 $f_{p}$ 函式來計算圓心和曲率。它用於計算幾乎所有沒有外圓的圓和硬圓。

函式 $f_{pm}$ 使用 $f_{p}$ 函式來計算曲率和 $f_{m}$ 函式來計算圓心。它用於計算硬圓。

函式 $f_{mm}$ 使用 $f_{m}$ 函式來計算圓心和曲率。它用於計算外圓。

用於填充間隙的函式

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

構造

初始階段：3 個相互相切的圓 = $(k_{ia},k_{ib},k_{ic})$
零階段：2 個與前一個圓相切的圓：外圓 $ck_{0out}$ 和內圓 $ck_{0in}$ ，其中：

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

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

第一階段：在之前的圓之間有 6 個空隙（= 曲線三角形）。三個空隙在圓周上（靠近外圓），三個空隙在中間（靠近內圓）。在每個空隙中放一個圓。

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

每個放在空隙中的新圓都會形成三個新的空隙（三叉）。

第二階段：有 18（6*3）個新的空隙。使用 $f_{pp}$ 函式在每個空隙中放置一個圓，除了一個硬圓。

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

後續階段：在每個階段的每個空隙中放置一個圓。在每個階段，只有一個硬圓，它在 $1c$ 空隙中。

程式設計師演算法

由於所有圓都內切於外圓，因此從外圓開始更容易。

選擇二次影像的大小 $side$
計算外圓的半徑

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

將外圓放置在笛卡爾平面的第一象限

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

計算 3 個相等的初始圓的半徑^[10]

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

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

找到初始圓的圓心。例如，將一個圓 $ia$ 放在上排，兩個圓 ( $ib$ 和 $ic$ ) 放在下排。

然後這些圓的圓心是等邊三角形的頂點，三角形邊長 $=2*r_{i}$ 因此

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

計算內圓 $0in$

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

"現在我們有一個可以開始的配置"^[11] 有 6 個空隙要填補。5 個空隙很容易填補，而一個空隙 $1c$ 很難填補。

圓之間的關係

{\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})

圓圈數量

第 n 階段新圓圈數量 = $2*3^{n}$
n 階段後圓圈總數 = $2+3^{n+1}$

其中 $n>=0$

例如

階段 -1 = 初始配置（3 個內圓）
階段 0 給出 2 個新圓（一個外圓和一個最內圓）。所有圓圈 = 5
階段 1 給出 6 個新圓（所有圓圈 = 11）
階段 2 給出 18 個新圓（所有圓圈 = 29）
階段 3 給出 54 個新圓（所有圓圈 = 83）
階段 4 給出 162 個新圓（所有圓圈 = 245）
階段 5 給出 486 個新圓（所有圓圈 = 731）
階段 6 給出 1 458 個新圓（所有圓圈 = 2 189）
階段 7 給出 4 374 個新圓（所有圓圈 = 6 563）
階段 8 給出 13 122 個新圓（所有圓圈 = 19 685）
階段 9 給出 39 366 個新圓（所有圓圈 = 59 051）
階段 10 給出 118 098 個新圓（所有圓圈 = 177 149）
階段 11 給出 354 294 個新圓（所有圓圈 = 531 443）
階段 12 給出 1 062 882 個新圓（所有圓圈 = 1 594 325）
階段 13 給出 3 188 646 個新圓（所有圓圈 = 4 782 971）
階段 14 給出 9 565 938 個新圓（所有圓圈 = 14 348 909）
階段 15 給出 28 697 814 個新圓（所有圓圈 = 43 046 723）
階段 16 給出 86 093 442 個新圓（所有圓圈 = 129 140 165）
階段 17 給出 258 280 326 個新圓（所有圓圈 = 387 420 491）
階段 18 給出 774 840 978 個新圓（所有圓圈 = 1 162 261 469）
階段 19 給出 2 324 522 934 個新圓（所有圓圈 = 3 486 784 403）
階段 20 給出 6 973 568 802 個新圓（所有圓圈 = 10 460 353 205）

可以使用 Maxima CAS 程式碼計算

 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), " )");

圓圈數量快速增長，導致大檔案尺寸問題。可以

僅使用最多 5-7 階段。這似乎足以在短時間內生成良好的影像。^[12]
不要繪製太小的圓圈（半徑小於畫素尺寸）^[13]
將 svg 檔案轉換為 png（例如使用 Image Magic : "convert a.svg a.png"）。然後，階段 12 svg 影像為 128 MB，而 png 檔案僅為 57 KB。
高興的是你得到了如此詳細的影像，並等待很長時間載入檔案，或者看看系統如何掛起（:-））

曲率

注意，只有一個曲率是用 $f_{\color {Red}m}$ 函式計算的。

它是 $k_{0out}$ 。它也是唯一的負曲率。

所有其他曲率都是用 $f_{p}$ 函式計算的。

階段 $n\,$	曲率
-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)\,$

示例程式

Haskell 和 EDSL Diagrams

影像和程式碼來自 Brent Yorgey^[14]
diagrams 程式碼^[15]

Java

Ken Stephenson 編寫的 CirclePack 程式

Maxima CAS

程式在上面解釋過，原始碼在 Apollonian_rc_6.svg 的描述中。

Common Lisp

要執行，請將 apollonian.lisp 檔案放到您的主目錄和 Emacs/slime 中。要進行下一步更改，請更改 *max-step*

載入檔案以執行程式碼（批處理模式）

(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

這是一個增量演算法，它維護著一個包含三個相互接觸的圓的佇列。主迴圈從佇列的開頭開始，移除一個三元組，用所需的圓形填充它，將生成的圓形新增到一個用於繪製為 SVG 的圓形列表中，並將新增圓形後產生的三個三元組新增到計算佇列的末尾。作為額外功能，這些圓形被著色，因此沒有兩個相鄰的圓形具有相同的顏色。

用法：將兩個 SVG 模板部分和中間的 Javascript 部分合併到一個 SVG 檔案中，並載入到網頁瀏覽器（或另一個支援 Javascript 或 ECMAScript 的 SVG 瀏覽器）中。

<?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()