Problem 199: Iterative Circle Packing

Project Euler Problem 199 Solution

EulerSolve provides an optimized solution for Project Euler Problem 199, Iterative Circle Packing, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary A unit circle contains three equal circles, each tangent to the other two and also tangent to the outer circle. The curved triangular gaps that remain are then filled recursively: in every gap bounded by three mutually tangent circles, insert the unique new circle tangent to all three boundary circles, and repeat this for a fixed number of generations. The problem asks for the fraction of the unit disk that is still uncovered after 10 iterations. The geometry looks two-dimensional, but the implementations never need explicit coordinates. The entire computation can be driven by curvatures alone, which turns the packing into a clean recursive area sum. Mathematical Approach The key idea is to represent every circle by its signed curvature $k=1/r$. Inner circles have positive curvature, while the enclosing unit circle is treated as a circle of curvature $-1$. With that convention, every local configuration is governed by Descartes' circle theorem....

Detailed mathematical approach

Problem Summary

Three circles of equal radius lie inside a larger circle so that each pair of the four circles is tangent and the inner circles do not overlap. Four curved gaps remain uncovered. In each iteration a circle of the largest possible size is placed in every gap, creating new gaps for the next iteration. After $3$ iterations there are $108$ gaps, and the fraction of the area of the large circle not covered by circles is $0.06790342$, rounded to eight decimal places. We must find this fraction after $10$ iterations, rounded to eight decimal places in the format x.xxxxxxxx. ^[1]

The whole solution in one sentence: by Descartes' circle theorem, the circle placed in a gap whose sides have the curvatures $a,b,c$, one over their radii and negative for the enclosing outer circle, has the curvature $a+b+c+2\sqrt{ab+bc+ca}$ and splits the gap into three new gaps, so a short recursion adds up the areas of all circles placed in ten iterations.

Mathematical Approach

We find the starting circles, count the gaps, identify the largest circle in a gap and compute it with Descartes' theorem.

1. Curvature instead of radius

Only a ratio of areas is asked for, so the outer circle gets the radius $1$ and the centre $O$. We describe every circle by its curvature, one over its radius; a circle of curvature $x$ covers the fraction $1/x^2$ of the disc. The outer circle, which encloses all others, gets the curvature $-1$, as Descartes' theorem in Section 5 requires. Each circle lies in its own gap, so no two overlap, and after $n$ iterations the uncovered fraction is $U_n=1-\sum 1/x^2$, summed over all circles inside.

2. The starting circles

Let $r$ be the radius of the three starting circles. They touch one another, so their centres form an equilateral triangle with side $2r$, centred at $O$ by symmetry. The distance from the centre of an equilateral triangle to a vertex is the side divided by $\sqrt3$, and each starting circle touches the outer circle, so $1-r=2r/\sqrt3$. Hence

$$\begin{aligned}r&=\frac{1}{1+2/\sqrt3}=2\sqrt3-3,\\k&=\frac{1}{r}=1+\frac{2}{\sqrt3},\end{aligned}$$

so $r\approx0.4641$, and $k\approx2.1547$ is their curvature; the same ratio $2\sqrt3-3$ between the sizes of the outer and the inner curvature appears for symmetric gaskets in the literature. ^[6] The starting circles cover $3r^2$ of the disc, so $U_0=1-3r^2$${}=36\sqrt3-62\approx0.35382907$.

3. How the gaps multiply

The four starting gaps are the central gap between the three starting circles and three outer gaps, each between two starting circles and the outer circle. A gap is a curved triangle with arcs of three circles that touch in pairs as its sides. The circle placed in it touches all three sides, as Section 4 shows, and its three points of contact cut the gap into the new disc and three smaller gaps, each bounded by two old sides and the new circle. Each iteration therefore triples the number of gaps: after $n$ iterations there are $4\cdot3^n$ gaps, $108$ for $n=3$ as in the statement, and $3+4\,(1+3+\dots+3^{n-1})$${}=2\cdot3^n+1$ circles.

4. The largest circle in a gap

Inside a gap, the largest circle centred at $P$ has as its radius the distance from $P$ to the nearest side: $|PC|-\rho$ for an inner circle with centre $C$ and radius $\rho$, and $1-|PO|$ for the outer circle. This distance is continuous on the closed and bounded gap, so it has a maximum: a largest circle exists. ^[5] It touches all three sides. For each side let $e$ be the unit vector along which the distance to it grows fastest, from $C$ towards $P$, or from $P$ towards $O$ for the outer circle; a small step of $P$ at an acute angle to $e$ increases that distance. So a circle touching at most one side can be moved off all sides and enlarged, and so can one touching two sides, by a step along $e_1+e_2$, unless $e_1=-e_2$. For two inner circles that puts $P$ on the segment between their centres, which the two touching circles cover. For the outer circle and an inner circle of radius $\rho$, centred at distance $1-\rho$ from $O$, it puts $P$ on the line through that centre, on the far side of $O$; at distance $t$ from $O$ the two distances $1-t$ and $t+1-2\rho$ are equal only for $t=\rho$, which gives a circle of radius $1-\rho$. But every later circle lies in a starting gap: the central one has the area $A_0=(\sqrt3-\pi/2)\,r^2\approx0.035$, the triangle of the centres minus three sectors of $60^\circ$, and each outer one $(\pi U_0-A_0)/3\approx0.359$, both below $\pi r^2\approx0.677$. So $\rho\le r$, and that circle would have an area of at least $\pi(1-r)^2\approx0.902$, more than any gap.

Three circles that touch one another at three different points are touched by exactly two further circles, their Soddy circles, one of which may be a straight line. ^[3] ^[4] One of them is known for every gap, its partner outside it: the outer circle for the central gap, the third starting circle for an outer gap, and, when a circle is placed in a gap with sides $a,b,c$, the side $c$ for the new gap with the sides $a$, $b$ and the new circle, and so on. The largest circle in the gap is the other one.

5. Descartes' circle theorem

Four circles that touch one another in pairs at six different points, with curvatures $k_1,k_2,k_3,k_4$, satisfy Descartes' circle theorem ^[2]

$$\begin{aligned}&(k_1+k_2+k_3+k_4)^2\\&\quad=2\,(k_1^2+k_2^2+k_3^2+k_4^2),\end{aligned}$$

where a circle that encloses the other three gets a negative curvature. The starting circles are a check: with $k_1=-1$ and $k_2=k_3=k_4=k$ the theorem reads $(3k-1)^2=2(1+3k^2)$, that is $3k^2-6k-1=0$, whose positive root is $k=1+2/\sqrt3$. For a gap with sides $a,b,c$ put $s=a+b+c$ and $q=ab+bc+ca$. Since $a^2+b^2+c^2=s^2-2q$, a fourth circle of curvature $x$ that touches all three satisfies

$$\begin{aligned}x^2-2sx+s^2-4q&=0,\\x&=s\pm2\sqrt{q}.\end{aligned}$$

The two roots belong to the new circle and the partner, and, as the coefficient $-2s$ shows, they add up to $2s$.

6. The plus sign

With $p$ the curvature of the partner, the new circle therefore has the curvature $x=2s-p$. The programs use the square root instead, so we show $p\lt s$, which makes $p$ the smaller root. For the central gap $p=-1\lt3k$, and for an outer gap $p=k\lt2k-1$, as $k\gt1$. If $p\lt s$ in some gap, then $x=s+(s-p)\gt s\gt c$, because $a+b\gt0$: at most one side is the outer circle, and an inner circle has a curvature above $1$. So the partner $c$ of the new gap with sides $a,b,x$ satisfies $c\lt x\lt a+b+x$, and likewise for the other two new gaps. By induction $p\lt s$ in every gap, and

$$\begin{aligned}x&=s+2\sqrt{q}\\&=a+b+c+2\sqrt{ab+bc+ca}.\end{aligned}$$

For three circles that touch from outside, this is the known rule that the plus sign gives the inner Soddy circle. ^[3]

7. Worked example: the first iteration

In the central gap $s=3k$ and $q=3k^2$, so $x=(3+2\sqrt3)\,k$${}=7+4\sqrt3\approx13.928$, and the other root, $(3-2\sqrt3)\,k=-1$, is the outer circle. By symmetry the new circle is centred at $O$ and touches circles whose centres lie at distance $1-r$, so its radius is $1-2r=7-4\sqrt3$, and indeed $(7+4\sqrt3)(7-4\sqrt3)=1$. In an outer gap $s=2k-1$ and $q=k^2-2k=1/3$, so $x=2k-1+2/\sqrt3$${}=1+2\sqrt3\approx4.464$, and the other root is $k$, the third starting circle. Geometrically, the new circle has its centre on the symmetry line of the gap, at a distance $t$ from $O$, and the radius $1-t$. A neighbouring starting centre lies at distance $1-r$ from $O$, at an angle of $60^\circ$ to that line, so by the law of cosines the two circles touch when

$$\begin{aligned}&t^2+(1-r)^2-t\,(1-r)\\&\quad=(1+r-t)^2.\end{aligned}$$

Hence $t=4r/(1+3r)$, and the radius is $(1-r)/(1+3r)=(2\sqrt3-1)/11$, which equals $1/(1+2\sqrt3)$. The first iteration covers $1/(7+4\sqrt3)^2$${}=97-56\sqrt3\approx0.00515$ and $3/(1+2\sqrt3)^2$${}=3(13-4\sqrt3)/121$${}\approx0.15054$, so

$$\begin{aligned}U_1&=U_0-(97-56\sqrt3)\\&\quad-\frac{3\,(13-4\sqrt3)}{121}\\&=\frac{11144\sqrt3-19278}{121}\\&\approx0.19813388.\end{aligned}$$

8. The recursion over the gaps

Let $F_d(a,b,c)$ be the fraction of the disc covered by the circles that $d$ iterations place in a gap with sides $a,b,c$. Then $F_0=0$ and, for $d\ge1$, with $x$ from Section 6,

$$\begin{aligned}F_d(a,b,c)&=\frac{1}{x^2}+F_{d-1}(a,b,x)\\&\quad+F_{d-1}(a,c,x)\\&\quad+F_{d-1}(b,c,x),\end{aligned}$$

since the circle $x$ comes first and the other $d-1$ iterations fill its three new gaps. By induction, $F_d$ depends only on the three curvatures, so the three outer gaps, all with the sides $-1,k,k$, receive the same area, and

$$\begin{aligned}U_n&=1-\frac{3}{k^2}-F_n(k,k,k)\\&\quad-3F_n(-1,k,k).\end{aligned}$$

For $n=0,1,\dots,5$ this gives, to eight decimals, $0.35382907$, $0.19813388$, $0.11413779$, $0.06790342$ as in the statement, $0.04167734$ and $0.02634471$.

9. Exact curvatures and slow corners

The rule $x=2(a+b+c)-p$ needs no square roots: from $-1$ and $k$ it builds every curvature as $u+vk$ with integers $u,v$, such as $7+4\sqrt3=1+6k$ and $1+2\sqrt3=3k-2$; the verification uses this exact form. It also shows why corners fill slowly. Each iteration puts one circle into the corner where two touching circles of curvatures $a$ and $b$ meet; it touches the previous circle of this chain and has the one before that as its partner, so

$$c_{j+1}=2\,(a+b+c_j)-c_{j-1},$$

and the second differences $c_{j+1}-2c_j+c_{j-1}=2(a+b)$ are constant, so the curvatures grow only quadratically. In the corner between the outer circle and a starting circle, $a=-1$, $b=k$ and $c_{-1}=c_0=k$, the other two starting circles, so $c_j=k+(k-1)\,j\,(j+1)$: $4.464$, $9.083$, $16.011$, $25.249$ and $36.796$ in iterations $1$ to $5$, the largest circles of these iterations. Only in the limit does the uncovered area vanish: the points never inside a circle form the Apollonian gasket, a set of area zero and fractal dimension about $1.3057$. ^[6]

10. Is double precision enough?

In double precision, with rounding to nearest, each operation on numbers of normal size has a relative error of at most $2^{-53}\approx1.1\cdot10^{-16}$, ^[9] but such errors can grow, for instance through cancellation in $q$ when $a=-1$. So we compare with exact values. Interval arithmetic carries guaranteed lower and upper bounds through every operation; ^[10] with the exact curvatures $u+vk$ of Section 9 it encloses each $U_n$, $n\le10$, in an interval narrower than $10^{-46}$. The double-precision recursion of the programs stays within $2.4\cdot10^{-16}$ of these values, and the value after $10$ iterations lies more than $4\cdot10^{-9}$ from the nearest point where its rounding to eight decimals would change, so the printed digits are certainly right.

How the Code Works

The C++ program ^[7] keeps depth, by default $10$, and the flag run_checkpoints in the struct Options. parse_arguments accepts --skip-checkpoints and --depth= followed by decimal digits, which parse_int_after_prefix collects into an int without an overflow check, so a huge number causes signed overflow, undefined behaviour in C++; ^[8] any other argument ends the program with exit code $1$. Unless skipped, run_checkpoints compares the results for depths $0$ and $1$ with the constants 0.3538290724795825 and 0.19813388055841763, that is $U_0$ and $U_1$, within $10^{-15}$; a failure gives exit code $2$.

fill_gap(k1, k2, k3, depth) returns $\pi F_d$ of Section 8, the covered area itself: at depth $0$ it returns $0$; otherwise it computes k4 by the formula of Section 6 and returns $\pi/k_4^2$ plus the calls for the gaps (k1, k2, k4), (k1, k3, k4) and (k2, k3, k4). uncovered_ratio takes $\pi$ as std::acos(-1.0L), $k=1+2/\sqrt3$ and $r=1/k$, adds to the area $3\pi r^2$ of the starting circles the call fill_gap(k, k, k, depth) and three times fill_gap(-1.0L, k, k, depth), and returns one minus the total divided by $\pi$, which main prints as a double with eight decimals.

The C++ code computes in long double, the 80-bit format with a 64-bit significand on many x86 compilers, but the 53-bit double format on the Apple arm64 machine used here. ^[8] Python nests the same recursion in solve, and Java has it in the static method fillGap; both use IEEE doubles, ^[9] fix the depth at $10$, run no checkpoints and print eight decimals, but Java's printf follows the default locale and writes a decimal comma under a German locale, for example. At depth $10$ the recursion is at most eleven calls deep.

Complexity and Verification

A call with depth $d\ge1$ makes three calls with depth $d-1$, so one top-level call makes $(3^{d+1}-1)/2$ calls and computes $(3^d-1)/2$ circles. With the factor $3$, two top-level calls suffice: for $d=10$, $177146$ calls and $59048$ computed circles, half of the $118096$ circles placed in the four gaps; the C++ checkpoints add $10$ calls. The time grows like $3^d$ and the memory like $d$; the C++ program finishes in a few milliseconds, Python and Java in a few hundredths of a second, mostly start-up time.

For this article the C++ program, built with GCC and with Clang, and the Python and Java programs were run: they print the same eight decimals and, on the test machine, compute the same double. Independent checks, not part of the published solutions:

The rule $x=2(a+b+c)-p$ in exact integer pairs $(u,v)$, with the interval arithmetic of Section 10, reproduces $U_0=36\sqrt3-62$, $U_1$, the $108$ gaps and the statement's $0.06790342$; all $118096$ gaps satisfy $p\lt s$ and $q\ge1/3$.
The same rule applied to curvature times centre, which the complex form of Descartes' theorem allows, ^[2] gives all $1459$ circles of six iterations with centres: each new circle touches its three sides to within $7\cdot10^{-16}$, no two circles overlap, and the midpoints of a $6000\times6000$ grid of squares give the uncovered fraction $0.067897$ after three iterations.
In the $26$ gaps filled in the first three iterations of the central gap and one outer gap, a grid search over $801\times801$ centres found no circle larger than the one touching all three sides.

Footnotes and References

The references below supply the background facts used above; the proofs of Sections 4 and 6, the worked example and all numbers are derived or computed in this article.

Project Euler 199 — Iterative Circle Packing. The official statement fills the four gaps between three equal, mutually tangent circles inside a larger circle with maximally sized circles, gives 108 gaps and the uncovered fraction 0.06790342 after 3 iterations, and asks for that fraction after 10 iterations, to eight decimals.
Wikipedia — Descartes' theorem. Four circles tangent to each other at six distinct points satisfy (k₁ + k₂ + k₃ + k₄)² = 2(k₁² + k₂² + k₃² + k₄²), with negative curvature for a circumscribing circle; k₄ = k₁ + k₂ + k₃ ± 2√(k₁k₂ + k₂k₃ + k₃k₁) gives the two circles tangent to a triple, and a complex form gives the centres.
MathWorld — Soddy Circles. For three pairwise externally tangent circles the plus branch gives the inner and the minus branch the outer Soddy circle, which can be a line.
Wikipedia — Problem of Apollonius. For three mutually tangent circles the problem has five solutions: the three given circles themselves and two further circles, Soddy's circles.
Wikipedia — Extreme value theorem. A continuous real function on a nonempty compact set, such as a closed bounded set in the plane, attains a maximum.
Wikipedia — Apollonian gasket. With three equal inner circles, the bounding circle's curvature is, in absolute value, 2√3 − 3 times theirs; the Hausdorff dimension is about 1.3056867. MathWorld — Apollonian Gasket: iterated inner Soddy circles leave a set of measure zero uncovered.
C++ — Euler199.cpp. The linked immutable revision, with fill_gap, uncovered_ratio and run_checkpoints, is the implementation record for this article.
cppreference — Fundamental types: long double need not be an IEEE-754 type; many x86 and x86-64 implementations use the x87 80-bit format, MSVC that of double. Arithmetic operators: signed integer overflow is undefined behaviour.
Wikipedia — Double-precision floating-point format: 53-bit significand, relative rounding error at most 2⁻⁵³ ≈ 1.11 × 10⁻¹⁶. Python maps its floats to binary64 on almost all platforms; the Java Language Specification, §4.2.3, associates double with binary64, and PrintStream.printf formats with the default locale.
Wikipedia — Interval arithmetic. Each value is a range with lower and upper bounds that also track rounding errors, so the computed ranges are guaranteed to contain the exact results.

Problem 199 source code

C++

#include <cmath>
#include <iomanip>
#include <iostream>
#include <string>

namespace {

struct Options {
    int depth = 10;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--depth=", options.depth)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.depth >= 0;
}

long double fill_gap(const long double k1,
                     const long double k2,
                     const long double k3,
                     const int depth) {
    if (depth == 0) {
        return 0.0L;
    }

    const long double k4 =
        k1 + k2 + k3 + 2.0L * std::sqrt(k1 * k2 + k2 * k3 + k3 * k1);
    const long double area = std::acos(-1.0L) / (k4 * k4);

    return area +
           fill_gap(k1, k2, k4, depth - 1) +
           fill_gap(k1, k3, k4, depth - 1) +
           fill_gap(k2, k3, k4, depth - 1);
}

long double uncovered_ratio(const int depth) {
    const long double pi = std::acos(-1.0L);
    const long double k = 1.0L + 2.0L / std::sqrt(3.0L);
    const long double r = 1.0L / k;

    long double covered = 3.0L * pi * r * r;
    covered += fill_gap(k, k, k, depth);
    covered += 3.0L * fill_gap(-1.0L, k, k, depth);

    return 1.0L - covered / pi;
}

bool run_checkpoints() {
    if (std::abs(uncovered_ratio(0) - 0.3538290724795825L) > 1e-15L) {
        std::cerr << "Checkpoint failed for depth=0" << '\n';
        return false;
    }
    if (std::abs(uncovered_ratio(1) - 0.19813388055841763L) > 1e-15L) {
        std::cerr << "Checkpoint failed for depth=1" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << std::fixed << std::setprecision(8)
              << static_cast<double>(uncovered_ratio(options.depth)) << '\n';
    return 0;
}

Python

# Problem 199: Iterative Circle Packing
# Apollonius gasket: 3 mutually tangent circles inside a larger circle.
# Find fraction of area not covered after 10 iterations.

import math

def solve():
    depth = 10

    def fill_gap(k1, k2, k3, d):
        if d == 0: return 0.0
        k4 = k1 + k2 + k3 + 2.0 * math.sqrt(k1*k2 + k2*k3 + k3*k1)
        area = math.pi / (k4*k4)
        return area + fill_gap(k1, k2, k4, d-1) + fill_gap(k1, k3, k4, d-1) + fill_gap(k2, k3, k4, d-1)

    pi = math.pi
    k = 1.0 + 2.0/math.sqrt(3.0)
    r = 1.0/k
    covered = 3.0 * pi * r * r
    covered += fill_gap(k, k, k, depth)
    covered += 3.0 * fill_gap(-1.0, k, k, depth)
    ratio = 1.0 - covered/pi
    print(f"{ratio:.8f}")

solve()

Java

public class Euler199 {
    static double fillGap(double k1, double k2, double k3, int depth) {
        if (depth == 0)
            return 0;
        double k4 = k1 + k2 + k3 + 2 * Math.sqrt(k1 * k2 + k2 * k3 + k3 * k1);
        double area = Math.PI / (k4 * k4);
        return area + fillGap(k1, k2, k4, depth - 1) + fillGap(k1, k3, k4, depth - 1) + fillGap(k2, k3, k4, depth - 1);
    }

    public static void main(String[] args) {
        int depth = 10;
        double k = 1 + 2 / Math.sqrt(3);
        double r = 1 / k;
        double covered = 3 * Math.PI * r * r;
        covered += fillGap(k, k, k, depth);
        covered += 3 * fillGap(-1, k, k, depth);
        System.out.printf("%.8f%n", 1 - covered / Math.PI);
    }
}