Problem 894: Spiral of Circles

View on Project Euler

Project Euler Problem 894 Solution

EulerSolve provides an optimized solution for Project Euler Problem 894, Spiral of Circles, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The problem asks for the area generated by an infinite self-similar spiral of tangent circles. The geometry is modeled in the complex plane by choosing a scale factor \(0<s<1\), a rotation angle \(t\), and a spacing factor \(d>0\), then placing circle \(i\) at center \(z_i=dq^i\) with radius \(r_i=s^i\), where \(q=se^{it}\). The job is to recover the correct geometric parameters from the tangency pattern and then sum the infinitely many curvilinear gaps. Mathematical Approach The implementation turns the spiral geometry into a small nonlinear system in two real variables. After solving that system numerically, it uses self-similarity to convert the infinite area into a geometric series. Step 1: Parametrize the Spiral Write $$q=se^{it},\qquad z_i=dq^i,\qquad r_i=s^i.$$ Because \(|q|=s\), moving from one circle to the next multiplies all lengths by \(s\) and rotates the picture by angle \(t\). This makes the packing self-similar: every deeper part of the spiral is a scaled copy of what came before. Step 2: Translate Tangency into Equations If circles \(i\) and \(i+m\) are tangent, then the distance between their centers equals the sum of their radii: $$|z_{i+m}-z_i|=r_i+r_{i+m}.$$ Substituting the model gives $$|d q^i(q^m-1)|=s^i+s^{i+m}.$$ After dividing by \(s^i\), the dependence on \(i\) disappears: $$d|1-q^m|=1+s^m.$$ This is the crucial simplification....

Detailed mathematical approach

Problem Summary

The problem asks for the area generated by an infinite self-similar spiral of tangent circles. The geometry is modeled in the complex plane by choosing a scale factor \(0<s<1\), a rotation angle \(t\), and a spacing factor \(d>0\), then placing circle \(i\) at center \(z_i=dq^i\) with radius \(r_i=s^i\), where \(q=se^{it}\). The job is to recover the correct geometric parameters from the tangency pattern and then sum the infinitely many curvilinear gaps.

Mathematical Approach

The implementation turns the spiral geometry into a small nonlinear system in two real variables. After solving that system numerically, it uses self-similarity to convert the infinite area into a geometric series.

Step 1: Parametrize the Spiral

Write

$$q=se^{it},\qquad z_i=dq^i,\qquad r_i=s^i.$$

Because \(|q|=s\), moving from one circle to the next multiplies all lengths by \(s\) and rotates the picture by angle \(t\). This makes the packing self-similar: every deeper part of the spiral is a scaled copy of what came before.

Step 2: Translate Tangency into Equations

If circles \(i\) and \(i+m\) are tangent, then the distance between their centers equals the sum of their radii:

$$|z_{i+m}-z_i|=r_i+r_{i+m}.$$

Substituting the model gives

$$|d q^i(q^m-1)|=s^i+s^{i+m}.$$

After dividing by \(s^i\), the dependence on \(i\) disappears:

$$d|1-q^m|=1+s^m.$$

This is the crucial simplification. For any fixed offset \(m\), the same tangency rule applies everywhere in the spiral.

Step 3: Use the Three Required Offsets

The intended configuration is pinned down by the tangent offsets \(m=1,7,8\). Therefore

$$d|1-q|=1+s,\qquad d|1-q^7|=1+s^7,\qquad d|1-q^8|=1+s^8.$$

Eliminating \(d\) with the first equation leaves two real equations in \(s\) and \(t\):

$$F_7(s,t)=(1+s)|1-q^7|-(1+s^7)|1-q|=0,$$

$$F_8(s,t)=(1+s)|1-q^8|-(1+s^8)|1-q|=0.$$

Since \(q^m=s^m e^{imt}\), the magnitude can be written explicitly as

$$|1-q^m|=\sqrt{1+s^{2m}-2s^m\cos(mt)}.$$

So the geometry has been reduced to a two-variable nonlinear system, which is exactly what Newton's method is designed to solve.

Step 4: Recover the Spiral Scale

Once \((s,t)\) is known, the remaining spacing factor follows immediately from the offset \(m=1\):

$$d=\frac{1+s}{|1-q|}.$$

The numerical solution reached by the implementations is approximately

$$s\approx 0.906331406148595,\qquad t\approx 0.826729539414059,\qquad d\approx 2.473989946674087.$$

At these values, the tangency identities for \(m=1,7,8\) hold to numerical precision, which confirms that the recovered spiral matches the intended pattern.

Step 5: Compute a Single Curvilinear Gap

For three mutually tangent circles with radii \(r_1,r_2,r_3\), the triangle formed by their centers has side lengths

$$a=r_2+r_3,\qquad b=r_1+r_3,\qquad c=r_1+r_2.$$

Its ordinary area is given by Heron's formula:

$$\Delta=\sqrt{p(p-a)(p-b)(p-c)},\qquad p=\frac{a+b+c}{2}.$$

If \(A_1,A_2,A_3\) are the corresponding angles of that center triangle, then the area trapped between the three circles is

$$\mathcal{C}(r_1,r_2,r_3)=\Delta-\frac{1}{2}\left(r_1^2A_1+r_2^2A_2+r_3^2A_3\right).$$

The angles come from the law of cosines, for example

$$A_1=\arccos\left(\frac{b^2+c^2-a^2}{2bc}\right),$$

and similarly for \(A_2\) and \(A_3\).

Step 6: Sum the Infinite Spiral by Self-Similarity

The implementations need two base curvilinear pieces:

$$A=\mathcal{C}(1,s,s^8),\qquad B=\mathcal{C}(1,s^7,s^8).$$

Every deeper layer of the spiral is similar to the previous one, with all lengths multiplied by \(s\). Areas therefore scale by \(s^2\), so the entire infinite sum is

$$A+B+s^2(A+B)+s^4(A+B)+\cdots=\frac{A+B}{1-s^2}.$$

This closed form is the reason the computation stays short even though the packing itself is infinite.

Worked Example

Using the converged value of \(s\), the two base pieces are approximately

$$A\approx 0.082310769808360,\qquad B\approx 0.055516558200733.$$

The common area ratio is

$$s^2\approx 0.8214366235.$$

Therefore

$$\frac{A+B}{1-s^2}\approx \frac{0.137827328009093}{0.1785633765}\approx 0.7718678168.$$

This matches the final numerical value returned by the implementations.

How the Code Works

The C++, Python, and Java implementations all follow the same structure. They start from a reasonable guess for \(s\) and \(t\), evaluate the two tangency equations, and approximate the \(2\times2\) Jacobian with small finite differences. Each Newton step solves that linear system, updates \((s,t)\), wraps the angle back into \([0,2\pi)\), and stops when the correction becomes negligible.

After the nonlinear solve, the implementation reconstructs \(d\) from the first tangency equation, evaluates the two base curvilinear areas using Heron's formula and sector subtraction, and returns the geometric-series value \((A+B)/(1-s^2)\). The C++ implementation also performs an explicit geometric sanity check: the required offsets \(1,7,8\) are tangent, while other nearby tested pairs do not overlap.

Complexity Analysis

The numerical problem has fixed size: two unknowns, two equations, two base area evaluations, and a bounded set of geometric checks. If \(I\) is the number of Newton iterations and \(V\) is the number of validated circle pairs, then the running time is \(O(I+V)\) and the memory usage is \(O(1)\). For the actual Project Euler instance, both \(I\) and \(V\) are small constants, so the method is effectively constant-time and constant-space.

Footnotes and References

  1. Problem page: Project Euler 894
  2. Newton's method: Wikipedia — Newton's method
  3. Circle packing: Wikipedia — Circle packing
  4. Heron's formula: Wikipedia — Heron's formula
  5. Law of cosines: Wikipedia — Law of cosines
  6. Geometric series: Wikipedia — Geometric series

Problem 894 source code

C++

#include <cassert>
#include <cmath>
#include <complex>
#include <iomanip>
#include <iostream>

using ld = long double;
using cd = std::complex<ld>;

static constexpr ld PI = 3.141592653589793238462643383279502884L;

static ld abs_one_minus_q_pow(ld s, ld t, int n) {
    ld rn = std::pow(s, static_cast<ld>(n));
    ld ang = static_cast<ld>(n) * t;
    cd qn = std::polar(rn, ang);
    return std::abs(cd(1, 0) - qn);
}

static std::pair<ld, ld> equations(ld s, ld t) {
    ld a1 = abs_one_minus_q_pow(s, t, 1);
    ld a7 = abs_one_minus_q_pow(s, t, 7);
    ld a8 = abs_one_minus_q_pow(s, t, 8);

    ld f7 = (1 + s) * a7 - (1 + std::pow(s, 7)) * a1;
    ld f8 = (1 + s) * a8 - (1 + std::pow(s, 8)) * a1;
    return {f7, f8};
}

static ld circular_triangle_area(ld r1, ld r2, ld r3) {
    ld a = r2 + r3;
    ld b = r1 + r3;
    ld c = r1 + r2;

    ld p = (a + b + c) / 2;
    ld tri = std::sqrt(std::max<ld>(0, p * (p - a) * (p - b) * (p - c)));

    ld A1 = std::acos((b * b + c * c - a * a) / (2 * b * c));
    ld A2 = std::acos((a * a + c * c - b * b) / (2 * a * c));
    ld A3 = std::acos((a * a + b * b - c * c) / (2 * a * b));

    return tri - 0.5L * (r1 * r1 * A1 + r2 * r2 * A2 + r3 * r3 * A3);
}

int main() {
    ld s = 0.9L;
    ld t = 0.83L;

    for (int it = 0; it < 80; ++it) {
        auto [f1, f2] = equations(s, t);

        ld h = 1e-12L;
        auto [f1s, f2s] = equations(s + h, t);
        auto [f1t, f2t] = equations(s, t + h);

        ld j11 = (f1s - f1) / h;
        ld j21 = (f2s - f2) / h;
        ld j12 = (f1t - f1) / h;
        ld j22 = (f2t - f2) / h;

        ld det = j11 * j22 - j12 * j21;
        assert(std::fabsl(det) > 1e-24L);

        ld ds = (-f1 * j22 + f2 * j12) / det;
        ld dt = (-j11 * f2 + j21 * f1) / det;

        s += ds;
        t += dt;

        while (t < 0) t += 2 * PI;
        while (t >= 2 * PI) t -= 2 * PI;

        if (std::fabsl(ds) + std::fabsl(dt) < 1e-20L) break;
    }

    auto [rf1, rf2] = equations(s, t);
    assert(std::fabsl(rf1) < 1e-12L);
    assert(std::fabsl(rf2) < 1e-12L);

    ld d = (1 + s) / abs_one_minus_q_pow(s, t, 1);
    cd q = std::polar(s, t);

    for (int i = 0; i < 20; ++i) {
        for (int j = i + 1; j < 40; ++j) {
            ld ri = std::pow(s, static_cast<ld>(i));
            ld rj = std::pow(s, static_cast<ld>(j));
            cd zi = d * std::pow(q, i);
            cd zj = d * std::pow(q, j);
            ld gap = std::abs(zi - zj) - (ri + rj);

            int diff = j - i;
            if (diff == 1 || diff == 7 || diff == 8) {
                assert(std::fabsl(gap) < 1e-10L);
            } else {
                assert(gap > -1e-10L);
            }
        }
    }

    ld A = circular_triangle_area(1, s, std::pow(s, 8));
    ld B = circular_triangle_area(1, std::pow(s, 7), std::pow(s, 8));
    ld total = (A + B) / (1 - s * s);

    std::cout << std::fixed << std::setprecision(10) << static_cast<double>(total) << '\n';
    return 0;
}

Python

import math
import cmath

PI = math.pi

def abs_one_minus_q_pow(s, t, n):
    rn = s ** n
    ang = n * t
    qn = cmath.rect(rn, ang)
    return abs(1.0 - qn)

def equations(s, t):
    a1 = abs_one_minus_q_pow(s, t, 1)
    a7 = abs_one_minus_q_pow(s, t, 7)
    a8 = abs_one_minus_q_pow(s, t, 8)

    f7 = (1.0 + s) * a7 - (1.0 + s ** 7) * a1
    f8 = (1.0 + s) * a8 - (1.0 + s ** 8) * a1
    return f7, f8

def circular_triangle_area(r1, r2, r3):
    a = r2 + r3
    b = r1 + r3
    c = r1 + r2

    p = (a + b + c) / 2.0
    val = p * (p - a) * (p - b) * (p - c)
    tri = math.sqrt(val) if val > 0 else 0.0

    A1 = math.acos((b * b + c * c - a * a) / (2.0 * b * c))
    A2 = math.acos((a * a + c * c - b * b) / (2.0 * a * c))
    A3 = math.acos((a * a + b * b - c * c) / (2.0 * a * b))

    return tri - 0.5 * (r1 * r1 * A1 + r2 * r2 * A2 + r3 * r3 * A3)

def solve():
    s = 0.9
    t = 0.83

    for _ in range(80):
        f1, f2 = equations(s, t)

        h = 1e-10
        f1s, f2s = equations(s + h, t)
        f1t, f2t = equations(s, t + h)

        j11 = (f1s - f1) / h
        j21 = (f2s - f2) / h
        j12 = (f1t - f1) / h
        j22 = (f2t - f2) / h

        det = j11 * j22 - j12 * j21
        
        ds = (-f1 * j22 + f2 * j12) / det
        dt = (-j11 * f2 + j21 * f1) / det

        s += ds
        t += dt

        while t < 0:
            t += 2 * PI
        while t >= 2 * PI:
            t -= 2 * PI

        if abs(ds) + abs(dt) < 1e-15:
            break

    A = circular_triangle_area(1.0, s, s ** 8)
    B = circular_triangle_area(1.0, s ** 7, s ** 8)
    total = (A + B) / (1.0 - s * s)

    return f"{total:.10f}"

if __name__ == "__main__":
    print(solve())

Java

public class Euler894 {

    static final double PI = Math.PI;

    static double absOneMinusQPow(double s, double t, int n) {
        double rn = Math.pow(s, n);
        double ang = n * t;
        double real = 1.0 - rn * Math.cos(ang);
        double imag = 0.0 - rn * Math.sin(ang);
        return Math.sqrt(real * real + imag * imag);
    }

    static class Pair {
        double f1, f2;

        Pair(double f1, double f2) {
            this.f1 = f1;
            this.f2 = f2;
        }
    }

    static Pair equations(double s, double t) {
        double a1 = absOneMinusQPow(s, t, 1);
        double a7 = absOneMinusQPow(s, t, 7);
        double a8 = absOneMinusQPow(s, t, 8);

        double f7 = (1.0 + s) * a7 - (1.0 + Math.pow(s, 7)) * a1;
        double f8 = (1.0 + s) * a8 - (1.0 + Math.pow(s, 8)) * a1;
        return new Pair(f7, f8);
    }

    static double circularTriangleArea(double r1, double r2, double r3) {
        double a = r2 + r3;
        double b = r1 + r3;
        double c = r1 + r2;

        double p = (a + b + c) / 2.0;
        double val = p * (p - a) * (p - b) * (p - c);
        double tri = val > 0 ? Math.sqrt(val) : 0.0;

        double A1 = Math.acos((b * b + c * c - a * a) / (2.0 * b * c));
        double A2 = Math.acos((a * a + c * c - b * b) / (2.0 * a * c));
        double A3 = Math.acos((a * a + b * b - c * c) / (2.0 * a * b));

        return tri - 0.5 * (r1 * r1 * A1 + r2 * r2 * A2 + r3 * r3 * A3);
    }

    public static String solve() {
        double s = 0.9;
        double t = 0.83;

        for (int it = 0; it < 80; ++it) {
            Pair f = equations(s, t);

            double h = 1e-10;
            Pair fs = equations(s + h, t);
            Pair ft = equations(s, t + h);

            double j11 = (fs.f1 - f.f1) / h;
            double j21 = (fs.f2 - f.f2) / h;
            double j12 = (ft.f1 - f.f1) / h;
            double j22 = (ft.f2 - f.f2) / h;

            double det = j11 * j22 - j12 * j21;

            double ds = (-f.f1 * j22 + f.f2 * j12) / det;
            double dt = (-j11 * f.f2 + j21 * f.f1) / det;

            s += ds;
            t += dt;

            while (t < 0)
                t += 2 * PI;
            while (t >= 2 * PI)
                t -= 2 * PI;

            if (Math.abs(ds) + Math.abs(dt) < 1e-15)
                break;
        }

        double A = circularTriangleArea(1.0, s, Math.pow(s, 8));
        double B = circularTriangleArea(1.0, Math.pow(s, 7), Math.pow(s, 8));
        double total = (A + B) / (1.0 - s * s);

        return String.format(java.util.Locale.US, "%.10f", total);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}