How the recursive structure of Apollonian gaskets can be described in order to be able to reproduce them?

Question

The classical Descartes-Soddy relationship between the signed curvatures $b_k$ ("b" for "bend") of 4 mutually tangent circles (Apollonian configuration):

$$\sum_{k=1}^4 b_k^2=\tfrac12 \left(\sum_{k=1}^4 b_k\right)^2\tag{1}$$

allows to obtain the curvature $b_4$, knowing $b_1,b_2,b_3$ by considering (1) as a quadratic equation in variable $b_4$. The fact that there are two solutions $b_4$ and $b'_4$ is in harmony with our experience. On fig. 1 below are represented 3 given mutually tangent circles in blue, and interior and exterior tangent circles to them in red).

Fig. 1.

From there, one obtains the radii $r_4=\tfrac{1}{|b_4|}$ and $r'_4=\tfrac{1}{|b'_4|}.$

But the centers $z_4$ and $z'_4$ (we work with complex notations) of the fourth circles are usually computed in a separate way.

In fact, the following unexpected formula (obtained only some 20 years ago) gives an easy way to obtain as well the centers of these circles:

$$\sum_{k=1}^4 (b_kz_k)^2=\tfrac12 \left(\sum_{k=1}^4 b_kz_k\right)^2\tag{2}$$

(please note the beautiful similarity with (1)!). One can find a proof of (2) with nD extensions by its discoverers in this well written document : https://arxiv.org/pdf/math/0101066.pdf

As before for relationship (1), $z_4$ is computed by considering (2) as a quadratic equation, this time with variable $z_4$ (or $b_4z_4$) giving the two centers $z_4$ and $z'_4$ (assuming that $b_4$ and $b'_4$ have been computed beforehand). See paragraph "complex Descartes theorem" in (https://en.wikipedia.org/wiki/Descartes%27_theorem).

Using formulas (1) et (2), I am able to make a certain number of steps by iterating initial step (figure 1); here is for example a second step where 6 new circles have been added to figure 1:

Fig. 2.

In the upsaid arxiv document, one finds in particular this figure :

Fig. 3 : A so-called "Apollonian gasket".

where the numbers figuring inside the disks are their unsigned curvatures.

I would like to "programmaticaly" reproduce this figure or similar figures but I am facing the difficulty to understand/manage its underlying recursive structure. Has somebody a hint ?

Valuable references : Indra's pearls (see chapter 7) : https://www.labri.fr/perso/mazoit/uploads/Book.pdf

https://www.americanscientist.org/article/a-tisket-a-tasket-an-apollonian-gasket

https://arxiv.org/ftp/arxiv/papers/0706/0706.0372.pdf

https://mathoverflow.net/q/88353

Other ones :

http://www.malinc.se/math/geometry/apolloniangasketen.php

https://arxiv.org/pdf/1309.3267.pdf

http://paulbourke.net/fractals/apollony/

And these two from the same site:

https://geometricolor.wordpress.com/2019/07/10/apollonian-gasket-as-a-spherical-fractal-with-tetrahedral-symmetry/

https://geometricolor.wordpress.com/2019/07/24/a-variant-of-the-apollonian-gasket-with-icosahedral-symmetry/

Remark 1 : The authors of the text mentionned at the beginning have published it in American Mathematical Monthly one year later: Jeffrey C. Lagarias, Colin L. Mallows, and Allan R. Wilks, Beyond the Descartes circle theorem, Amer. Math. Monthly 109 (2002), no. 4, 338–361.

Remark 2 : The term "Circles of Apollonius" may be misleading ; there are other circles having this name as explained here.

Remark 3: There are strong connections with Farey sequences (see in particular the last figure of this article).

Answer 1

I wrote a Mathematica function to generate such gaskets, and the way it navigates the net of circles is based on a function that returns the indices of the $3$ predecessors for the circle with index $n$:

Pred[n_Integer] := If[n < 6, {1, 2, 3}, Module[{q = Quotient[n, 3] + 2, p}, p = Pred[q]; Append[ Switch[Mod[n, 3], 0, p[[{1, 2}]], 1, p[[{1, 3}]], 2, p[[{2, 3}]]], q]]]

The circles $1$-$5$ are easy to generate explicitly from the first $3$ curvatures, $a\lt0$ and $b,c\gt0$:

$d=a+b+c-2\,\mathrm{Disc}(a,b,c)$
$e=a+b+c+2\,\mathrm{Disc}(a,b,c)$

$\text{circle }1:\left(\left(-\frac1a,0\right),\frac1a\right)$
$\text{circle }2:\left(\left(\frac1b,0\right),\frac1b\right)$
$\text{circle }3:\left(\left(\frac{b-a}{(a+b)c},\frac{2\,\mathrm{Disc}(a,b,c)}{(a+b)c}\right),\frac1c\right)$
$\text{circle }4:\left(\left(\frac{b-a}{(a+b)d},-\frac{2\,\mathrm{Disc}(a,b,d)}{(a+b)d}\right),\frac1d\right)$
$\text{circle }5:\left(\left(\frac{b-a}{(a+b)e},\frac{2\,\mathrm{Disc}(a,b,e)}{(a+b)e}\right),\frac1e\right)$

where $\mathrm{Disc}(a,b,c)=\sqrt{ab+bc+ca}$.

Then circles $\ge6$ can be computed using Pred[n] and a function that takes $3$ circles and returns the smaller circle touching all $3$:

NextCircle[a_Circle, b_Circle, c_Circle] := Module[{wa, ka = Curv[a], wb, kb = Curv[b], wc, kc = Curv[c], kd}, kd = ka + kb + kc + 2 Disc[ka, kb, kc]; wa = ka Disc[kb, kc, kd]; wb = kb Disc[kc, kd, ka]; wc = kc Disc[kd, ka, kb]; Circle[(wa Cent[a] + wb Cent[b] + wc Cent[c])/(wa + wb + wc), 1/kd]]

where Curv[c] returns the reciprocal of the radius of $c$ and Cent[c] returns the center of $c$.

Here is the result for $(a,b,c)=(-9,14,26)$:

Answer 2

Start off with the triple of circles given by the outer circle of curvature $1$ and the two circles with curvature $2$. In each step, find the interior tangent circle to a triple and generate three new triples, each containing the new circle and a pair of circles from the previous triple. To draw the circles in order of increasing curvature, immediately calculate the curvature for each triple as you generate it and keep them in a data structure sorted by curvature so that you can always process the smallest curvature remaining.

In the first step you'll have the two symmetric circles of curvature $3$ instead of an interior and an exterior one. You can save half the calculations by arbitrarily choosing one of these and then for each circle drawing both the circle and its vertical mirror image.