Apollonian Gaskets in Python

Wed 2014-08-27

Translations: de

An Apollonian Gasket with D3-Symmetry

Figure 1: An Apollonian Gasket with D3-Symmetry

Apollonian Gaskets are fractals that can be generated from three mutually tangent circles. From these, more circles which fill the enclosing circle can be calculated recursively.

I implemented this in python as a command line program that saves those as svg images. Below I will explain the math behind it and show some images.

An online version is also included so you can try it yourself!

As this turned out to be quite a lengthy article, here is a

1   How it works

The process of generating an Apollonian gasket roughly works like this:

We start with a triple of circles, each of which touches the other two from the outside. Now we try to find a fourth circle that is also tangent to each of the three. It's easy to see that there are two possibilities for this: Either it lies in the middle of the three given circles (externally tangent) or it encloses them (internally tangent).

Three mutually tangent circles, two alternatives for a fourth

Figure 2: Three and four mutually tangent circles

Left: Three circles with different radii. Each is tangent to the other two.

Right: Two possibilities for a fourth circle that is tangent to the first three. Externally tangent (pink) or internally tangent (green).

For now we want take the enclosing one.

There are two equations relating the curvatures and centers of four mutually tangent circles. That means, knowing our three circles from the start, we just plug in their values and get our enclosing circle.

Now comes the clever bit: The equations are quadratic polynomials, so we can use Vieta's formulas on them. This is very useful, because quadratic polynomials have two solutions (in this case: the inner and outer circle). With Vieta's formulas it is trivially easy to get the second solution if we already know the first.

Applied to our four circles that means:

  1. Considering our three circles form the start as "fixed" and the outer one as the first solution, Vieta gives us the small circle inside the three.
  2. Considering one of the starting circles as the first solution and keeping the rest fixed, Vieta gives another new circle. Because we started out with three circles, there are three candidates to plug into Vieta as the first solution, so we also get three new circles.
Four possibilities of applying Vieta's formulas

Figure 3: Four possibilities of applying Vieta's formulas.

The grey circles are fixed, the green circle is the first solution we already know about.

Pink is the second solution given by Vieta.

In total we have added four new circles.

Four new circles are added to the Apollonian gasket.

Figure 4: Four new circles are added to the Apollonian gasket.

This can be continued with our shiny new circles now added to the gasket: Find four mutually tangent ones, in turn consider each of them as the first of Vieta's solutions and get the second solution. Rinse and repeat.

New circles after second iteration of the Vieta process.

Figure 5: New circles after second iteration of the Vieta process.

When we think the gasket is full enough, we fill in some nice colors and save it to a file.

Now that you have an idea what this is about, you might want to read the rest of the article which details all the maths and programming aspects (which are not that bad, actually) or you may rather skip to the end of this article and try the Apollonian Gasket Generator for yourself: Try it!.

2   The Math

The math stuff can be broken into some smaller individual problems.

2.1   Why it works: Vieta revisited

If you are rusty on Vieta's formulas, here is a quick recap:

Consider the quadratic polynomial \(f(x)=x^2 + px + q\) with the solutions \(x_1\) and \(x_2\) of \(f(x)=0\). Vieta's formulas state that

\begin{equation} \label{eq:vietaplus} x_1 + x_2 = -p \end{equation}

\begin{equation} x_1\cdot x_2 = q \end{equation}

If we happen to know \(x_1\) already, \(x_2\) can be trivially computed as \(x_2=-p-x_1\).

Note that we do not use \(x_2=q/x_1\) as \(x_1\) might be \(0\) and the whole thing could explode in our faces with a "division by zero" error.

2.2   Representing Circles

There are (among others) two ways to represent a circle:

  1. Center and radius
  2. Center and curvature

The center is straightforward and can be thought of as a point in the 2D-plane or a complex number.

The radius as half the diameter is not even remotely difficult, only thing you might not already have heard is the curvature. Curvature is simply the reciprocal of the radius.

A circle with infinite radius has the curvature 0 and is, in fact, a straight line. On the other hand, infinite curvature means a radius of zero -- a point.

For us, it will be convenient to think of a circle as a complex number for the center and a curvature.

2.3   Descartes and his Disciples

Although the problem of the tangent circles was first described (and solved) by Apollonius of Perga (about 200 BC), doing this by hand with a compass and a straightedge is quite tedious still today.

The first step to an algorithmic solution was made by Descartes in 1643. He proved that the curvatures \(b_1, b_2, b_3, b_4\) of four mutually tangent circles satisfy

\begin{equation} \label{eq:descartes} b_1^2 + b_2^2 + b_3^2 + b_4^2 = \frac{1}{2}(b_1+b_2+b_3+b_4)^2 \end{equation}

This is known as Descartes' theorem.

Note that for this to work an enclosing (internally tangent) circle has to have a negative curvature.

Let's pretend we already know three circles' curvatures \(b_1,b_2,b_3\) and need to calculate \(b_4\). For this solve for \(b_4\), and after some rearranging of ugly long lines we get

\begin{equation} \label{eq:descartesb4} b_4 = b_1+b_2+b_3\pm 2\sqrt{b_1b_2 + b_2b_3 + b_3b_1} \end{equation}

But this only gives us the curvatures. Luckily, in 2001 Wilks et al discovered a very similar equation for the centers of Apollonian circles.

Let \(b_1,b_2,b_3,b_4\) be the curvatures and \(z_1,z_2,z_3,z_4\) the centers (as complex numbers) of four mutually tangent circles. They satisfy

\begin{equation} \label{eq:wilks} (b_1z_1)^2 + (b_2z_2)^2 + (b_3z_3)^2 + (b_4z_4)^2 = \frac{1}{2} (b_1z_1 + b_2z_2 + b_3z_3 + b_4z_4)^2 \end{equation}

Like Descartes' theorem, this can be solved for \(b_4z_4\), and with \(b_4\) from equation (\(\ref{eq:descartesb4}\)) it is a simple division to compute \(z_4\).

\begin{equation} \label{eq:wilksb4z4} b_4z_4 = b_1z_1+b_2z_2+b_3z_3\pm 2\sqrt{b_1b_2z_1z_2 + b_2b_3z_2z_3 + b_3b_1z_3z_1} \end{equation}

This does not look really elegant (square roots tend to evaluate to ugly irrational numbers, not to speak of square roots of complex numbers), but these equations are only needed one time: When we produce an enclosing circle for the three starting circles.

2.4   Vieta to the rescue!

Now that we have our four mutually tangent circles we can begin using one of Vieta's formulas (\(\ref{eq:vietaplus}\)) on them.

Let's first look at the curvatures:

Suppose we know values for \(b_i\) with \(i=1\dots4\) and are looking for the second solution for \(b_1\) that satisfies Descartes' theorem. We can rewrite it in a form more resembling the quadratic polynomial from (\(\ref{eq:vietaplus}\)).

First note that Descartes' therorem (\(\ref{eq:descartes}\)) can be written as \(2(b_1^2 + b_2^2 + b_3^2 + b_4^2) =(b_1+b_2+b_3+b_4)^2\), i.e. multiplied by two. From there we have:

\begin{eqnarray*} 2(b_1^2 + b_2^2 + b_3^2 + b_4^2) &=& (b_1+b_2+b_3+b_4)^2 \\ &=& b_1^2+b_2^2+b_3^2+b_4^2 +2b_1b_2 +2b_1b_3 +2b_1b_4 +2b_2b_3+2b_2b_4 +2b_3b_4 \\ &=& (b_1^2 + b_2^2 + b_3^2 + b_4^2) + 2b_1(b_2+b_3+b_4) +2(b_2b_3+b_2b_4 +b_3b_4) \end{eqnarray*}

Rearranging stuff to the left side we get

\begin{eqnarray*} b_1^2 + b_2^2 + b_3^2 + b_4^2 -2b_1(b_2+b_3+b_4)-2(b_2b_3+b_2b_4+b3_b4)&=& 0 \\ b_1^2 +\underbrace{(-2(b_2+b_3+b_4))}_{p}b_1 + \underbrace{(-2(b_2b_3+b_2b_4+b3_b4) + b_2^2 + b_3^2 + b_4^2)}_q&=&0 \end{eqnarray*}

Hooray, we're nearly there! Now we have something looking like a quadratic polynomial that we can feed to Vieta. Let \(b_{1,old},b_{1,new}\) be the solutions to \(b_1^2+pb_1 + q = 0\). Then we see from (\(\ref{eq:vietaplus}\)) that \(b_{1,old} + b_{1,new}=-p\).

\begin{equation} \label{eq:descartesvieta} b_{1,new}=-p - b_{1,old} = 2(b_2 + b_3 + b_4) - b_{1,old} \end{equation}

And that's it. With this equation, we can calculate the alternative fourth circle's curvature from the known four.

With the other equation, for the center points, it works the same, just replace \(b_i\) with \(b_iz_i\). In the end we get

\begin{equation} \label{eq:wilksvieta} z_{1,new}=\frac{2(b_2z_2 + b_3z_3 + b_4z_4) - b_{1,old}z_{1,old}}{b_{1,new}} \end{equation}

2.5   How to find starting triples

Now that we know how to proceed from four known circles to new ones we're only missing a method to come up with these three in the first place.

Finding appropriate circles is not that difficult either. Our goal is a method that, when given three curvatures \(b_1,b_2,b_3\), returns the curvature of a fourth, enclosing, circle and the centers of all of them.

  1. We will do geometry with the radii of the circles, so let \(r_i=1/b_i\).

  2. The first circle is easy. We stick it in the origin at \((0,0)\).

    First circle at the origin.
  3. No worries about the second circle either, just place it to the right of the first at \((r_1+r_2,0)\).

    Second circle to the right of number one.
  4. The third one is a bit more tricky. We want to place it below the first and second, but how do we get the coordinates of \((m_x,m_y)\)?

    Third circle below. What are the coordinates?

    To solve this, we take a look just at the triangle defined by the three radii.

    Triangle of the three radii.

    We can say that \((m_x,m_y)=(x,y)\) in the sketch above. We can further say that by Pythagoras' theorem the following must be the case:

    \begin{eqnarray*} y^2 + x^2&=& (r_1+r_3)^2 \\ y^2 + (r_1 + r_2 -x)^2 &=& (r_2+r_3)^2 \end{eqnarray*}

    Solving for \(x\) and \(y\) we get:

    \begin{eqnarray} \label{eq:thirdx} x = \frac{r_1^2 + r_1r_3 + r_1r_2 - r_2r_3}{r_1+r_2} \\ \label{eq:thirdy} y = \sqrt{(r_1+r_3)^2-x^2} \end{eqnarray}

  5. The fourth circle is easy again, we just feed everything about circles #1 to #3 into our equations (\(\ref{eq:descartesb4}\)) and (\(\ref{eq:wilksb4z4}\)) from earlier and get our last circle. Note that in place of the \(\pm\) we take a \(-\), because we want an enclosing circle, which has a negative curvature.

2.6   When it breaks

Note that it is not possible to generate Apollonian gaskets from all starting curvatures, for example, when the enclosing circle would need to have zero curvature to be tangent to the three.

One combination with this property is \(b_1=1,b_2=1,b_3=4\), or more generally any combination that satisfies \(b_3=2\sqrt{b_1\cdot b_2}+b_1+b_2\). In these cases equation (\(\ref{eq:descartesb4}\)) gives \(b_4=0\).

3   The Code

You can find the code on github.

First of all, we need to store the information about the circles. Unsurprisingly, the Circle class takes care of this.

 1 class Circle(object):
 2     """
 3     A circle represented by center point as complex number and radius.
 4     """
 5     def __init__ ( self, mx, my, r ):
 6         """
 7         @param mx: x center coordinate
 8         @param my: y center coordinate
 9         @param r: radius
10         """
11         self.r = r
12         self.m = (mx +my*1j)
14     def curvature (self):
15         return 1/self.r

Next, we want to get the enclosing circle for three mutually tangent circles. This is a straightforward implementation of equation (\(\ref{eq:descartesb4}\)) and (\(\ref{eq:wilksb4z4}\)).

 1 def outerTangentCircle( circle1, circle2, circle3 ):
 2     """
 3     Takes three externally tangent circles and calculates the fourth one enclosing them.
 4     """
 5     cur1 = circle1.curvature()
 6     cur2 = circle2.curvature()
 7     cur3 = circle3.curvature()
 8     m1 = circle1.m
 9     m2 = circle2.m
10     m3 = circle3.m
11     cur4 = -2 * sqrt( cur1*cur2 + cur2*cur3 + cur1 * cur3 ) + cur1 + cur2 + cur3
12     m4 = ( -2 * sqrt( cur1*m1*cur2*m2 + cur2*m2*cur3*m3 + cur1*m1*cur3*m3 ) + cur1*m1 + cur2*m2 + cur3*m3 ) /  cur4
13     circle4 = Circle( m4.real, m4.imag, 1/cur4 )
15     return circle4

We can directly use this if we comfortably want to generate the first four circles from three given radii. And while we are at it, we also implement equations (\(\ref{eq:thirdx}\)) and (\(\ref{eq:thirdy}\)) for the center of circle #3/4 (slightly other numbering in the code than in the equations).

 1 def tangentCirclesFromRadii( r2, r3, r4 ):
 2     """
 3     Takes three radii and calculates the corresponding externally
 4     tangent circles as well as a fourth one enclosing them. The
 5     enclosing circle is the first one.
 6     """
 8     circle2 = Circle( 0, 0, r2 )
 9     circle3 = Circle( r2 + r3, 0, r3 )
10     m4x = (r2*r2 + r2*r4 + r2*r3 - r3*r4) / (r2 + r3)
11     m4y = sqrt( (r2 + r4) * (r2 + r4) - m4x*m4x )
12     circle4 = Circle( m4x, m4y, r4 )
13     circle1 = outerTangentCircle( circle2, circle3, circle4 )
14     return ( circle1, circle2, circle3, circle4 )

The more interesting parts are next: Using Vieta to jump from circles we already know to new solutions:

 1 def secondSolution( fixed, c1, c2, c3 ):
 2     """
 3     If given four tangent circles, calculate the other one that is tangent
 4     to the last three.
 6     @param fixed: The fixed circle touches the other three, but not
 7     the one to be calculated.
 9     @param c1, c2, c3: Three circles to which the other tangent circle
10     is to be calculated.
11     """
13     curf = fixed.curvature()
14     cur1 = c1.curvature()
15     cur2 = c2.curvature()
16     cur3 = c3.curvature()
18     curn = 2 * (cur1 + cur2 + cur3) - curf
19     mn = (2 * (cur1*c1.m + cur2*c2.m + cur3*c3.m) - curf*fixed.m ) / curn
20     return Circle( mn.real, mn.imag, 1/curn )

The ApollonianGasket class holds the starting circles and a list of all generated circles.

1 class ApollonianGasket(object):
2     """
3     Container for an Apollonian Gasket.
4     """
5     def __init__(self, c1, c2, c3):
6         self.start = tangentCirclesFromRadii( 1/c1, 1/c2, 1/c3 )
7         self.genCircles = list(self.start)

And finally the most interesting part: Recursively generating the new circles from the starting set.

 1 def recurse(self, circles, depth, maxDepth):
 2     """Recursively calculate the smaller circles of the AG up to the
 3     given depth. Note that for depth n we get 2*3^{n+1} circles.
 5     @param maxDepth: Maximal depth of the recursion.
 6     @type maxDepth: int
 8     @param circles: 4-Tuple of circles for which the second
 9     solutions are calculated
10     @type circles: (L{Circle}, L{Circle}, L{Circle}, L{Circle})
12     @param depth: Current depth
13     @type depth: int
14     """
15     if( depth == maxDepth ):
16         return
17     (c1, c2, c3, c4) = circles
18     if( depth == 0 ):
19         # First recursive step, this is the only time we need to
20         # calculate 4 new circles.
21         del self.genCircles[4:]
22         cspecial = secondSolution( c1, c2, c3, c4 )
23         self.genCircles.append( cspecial )
24         self.recurse( (cspecial, c2, c3, c4), 1, maxDepth )
26     cn2 = secondSolution( c2, c1, c3, c4 )
27     self.genCircles.append( cn2 )
28     cn3 = secondSolution( c3, c1, c2, c4 )
29     self.genCircles.append( cn3 )
30     cn4 = secondSolution( c4, c1, c2, c3 )
31     self.genCircles.append( cn4 )
33     self.recurse( (cn2, c1, c3, c4), depth+1, maxDepth )
34     self.recurse( (cn3, c1, c2, c4), depth+1, maxDepth )
35     self.recurse( (cn4, c1, c2, c3), depth+1, maxDepth )

For each circle we construct the second solution by Vieta's formula and recurse with the original circle replaced by the new one. By sticking it in the first position of the circles quadruple and not recursing for the first circle (except in the first call of the function) we avoid that the second solution of the second solution is calculated, as this would again be the "grandparent" circle.

The remainder of the code just wraps all of this stuff in a command line program and adds some helper functions for svg output (ugly) and coloring.

4   Try it!

Another Apollonian gasket

Figure 6: Green is also a nice color.

In addition to the command line program I wrote a small online version. You can play around with different starting curvatures, recursion depths and colours.

You can find it at http://lsandig.org/cgi-bin/apollon/index.cgi


This text by Ludger Sandig is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.