Problem 727: Triangle of Circular Arcs

Project Euler Problem 727 Solution

EulerSolve provides an optimized solution for Project Euler Problem 727, Triangle of Circular Arcs, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Take three pairwise externally tangent circles with integer radii $a$, $b$, and $c$, where $$1 \le a \lt b \lt c \le 100,\qquad \gcd(a,b,c)=1.$$ Let $D$ be the circumcenter of the triangle formed by the three tangency points of the original circles, and let $E$ be the center of the inner Soddy circle tangent to all three of them. For each primitive triple define $$d(a,b,c)=\operatorname{dist}(D,E).$$ The goal is to average $d(a,b,c)$ over all primitive triples in the range. The implementation therefore needs an exact constant-time geometric formula for one triple and then a complete enumeration of all admissible triples. Mathematical Approach The geometry becomes straightforward once the three original circle centers are placed in coordinates. Every later quantity in the implementation is derived from that coordinate model. Step 1: Place the three circle centers Let $O_a$, $O_b$, and $O_c$ be the centers of the circles of radii $a$, $b$, and $c$. Because the circles are externally tangent, the distances between centers are $$|O_aO_b|=a+b,\qquad |O_aO_c|=a+c,\qquad |O_bO_c|=b+c.$$ Choose coordinates $$O_a=(0,0),\qquad O_b=(a+b,0).$$ The third center lies above the $x$-axis....

Detailed mathematical approach

Problem Summary

Take three pairwise externally tangent circles with integer radii $a$, $b$, and $c$, where

$$1 \le a \lt b \lt c \le 100,\qquad \gcd(a,b,c)=1.$$

Let $D$ be the circumcenter of the triangle formed by the three tangency points of the original circles, and let $E$ be the center of the inner Soddy circle tangent to all three of them. For each primitive triple define

$$d(a,b,c)=\operatorname{dist}(D,E).$$

The goal is to average $d(a,b,c)$ over all primitive triples in the range. The implementation therefore needs an exact constant-time geometric formula for one triple and then a complete enumeration of all admissible triples.

Mathematical Approach

The geometry becomes straightforward once the three original circle centers are placed in coordinates. Every later quantity in the implementation is derived from that coordinate model.

Step 1: Place the three circle centers

Let $O_a$, $O_b$, and $O_c$ be the centers of the circles of radii $a$, $b$, and $c$. Because the circles are externally tangent, the distances between centers are

Choose coordinates

$$O_a=(0,0),\qquad O_b=(a+b,0).$$

The third center lies above the $x$-axis. By projecting onto the base segment and applying the law of cosines, its coordinates are

$$x_c=\frac{(a+c)^2+(a+b)^2-(b+c)^2}{2(a+b)},$$

$$y_c=\sqrt{(a+c)^2-x_c^2},$$

$$O_c=(x_c,y_c).$$

Step 2: Build the tangency triangle and locate $D$

The tangency point between the circles of radii $a$ and $b$ lies on segment $O_aO_b$, at distance $a$ from $O_a$ and $b$ from $O_b$. Therefore

$$T_{ab}=O_a+\frac{a}{a+b}(O_b-O_a)=(a,0).$$

By the same ratio argument, the other two tangency points are

$$T_{ac}=O_a+\frac{a}{a+c}(O_c-O_a)=\frac{a}{a+c}O_c,$$

$$T_{bc}=O_b+\frac{b}{b+c}(O_c-O_b).$$

Point $D$ is the circumcenter of triangle $T_{ab}T_{ac}T_{bc}$. For any non-collinear points $T_i=(x_i,y_i)$, the circumcenter is given by

$$\Delta=2\bigl(x_1(y_2-y_3)+x_2(y_3-y_1)+x_3(y_1-y_2)\bigr),$$

$$D_x=\frac{(x_1^2+y_1^2)(y_2-y_3)+(x_2^2+y_2^2)(y_3-y_1)+(x_3^2+y_3^2)(y_1-y_2)}{\Delta},$$

$$D_y=\frac{(x_1^2+y_1^2)(x_3-x_2)+(x_2^2+y_2^2)(x_1-x_3)+(x_3^2+y_3^2)(x_2-x_1)}{\Delta}.$$

This is the exact formula used in the C++, Python, and Java implementations.

Step 3: Compute the inner Soddy circle radius

Let the curvatures of the three given circles be

$$k_1=\frac1a,\qquad k_2=\frac1b,\qquad k_3=\frac1c.$$

For the inner tangent circle, Descartes' circle theorem gives the positive fourth curvature

$$k_4=k_1+k_2+k_3+2\sqrt{k_1k_2+k_2k_3+k_3k_1}.$$

Hence the inner Soddy radius is

$$r=\frac1{k_4}.$$

If $E=(x_E,y_E)$ denotes the center of that circle, tangency implies

$$|EO_a|=a+r,\qquad |EO_b|=b+r,\qquad |EO_c|=c+r.$$

Step 4: Solve directly for the center $E$

Square the three distance conditions. From $O_a$ we get

$$x_E^2+y_E^2=(a+r)^2.$$

From $O_b=(a+b,0)$ we get

$$(x_E-(a+b))^2+y_E^2=(b+r)^2.$$

Subtracting these equations cancels the quadratic terms and leaves a linear equation:

$$2(a+b)x_E=(a+b)^2+(a+r)^2-(b+r)^2.$$

$$x_E=\frac{(a+b)^2+(a+r)^2-(b+r)^2}{2(a+b)}.$$

Now use the equation for $O_c=(x_c,y_c)$:

$$(x_E-x_c)^2+(y_E-y_c)^2=(c+r)^2.$$

Subtract the equation for $O_a$ again to obtain

$$2x_cx_E+2y_cy_E=x_c^2+y_c^2+(a+r)^2-(c+r)^2,$$

hence

$$y_E=\frac{x_c^2+y_c^2+(a+r)^2-(c+r)^2-2x_cx_E}{2y_c}.$$

The center of the inner Soddy circle is therefore obtained from two linear equations after one square root for $r$.

Step 5: Restrict to primitive triples and average

The outer search runs over all increasing triples

$$1 \le a \lt b \lt c \le 100$$

and keeps only those with

$$\gcd(a,b,c)=1.$$

This primitive condition removes scaled copies of the same similarity class. Indeed, if every radius is multiplied by a factor $t$, then every center, every tangency point, the circumcenter $D$, and the Soddy center $E$ all scale by $t$, so

$$d(ta,tb,tc)=t\,d(a,b,c).$$

After summing $d(a,b,c)$ over all primitive triples, the final result is just the arithmetic mean of those values.

Worked Example: $(a,b,c)=(1,2,3)$

This example is especially clean because the center triangle has side lengths $3$, $4$, and $5$. Therefore

$$O_a=(0,0),\qquad O_b=(3,0),\qquad O_c=(0,4).$$

The three tangency points are

$$T_{ab}=(1,0),\qquad T_{ac}=(0,1),\qquad T_{bc}=\left(\frac95,\frac85\right).$$

The circumcenter of triangle $T_{ab}T_{ac}T_{bc}$ is

$$D=(1,1).$$

The three curvatures are $1$, $1/2$, and $1/3$, so

$$k_1k_2+k_2k_3+k_3k_1=\frac12+\frac16+\frac13=1.$$

Descartes then yields

$$k_4=1+\frac12+\frac13+2=\frac{23}{6},\qquad r=\frac{6}{23}.$$

Substituting into the linear formulas gives

$$E=\left(\frac{21}{23},\frac{20}{23}\right).$$

Finally,

$$d(1,2,3)=|DE|=\sqrt{\left(1-\frac{21}{23}\right)^2+\left(1-\frac{20}{23}\right)^2}=\frac{\sqrt{13}}{23}\approx 0.1567630989,$$

which matches the sample value checked by the implementation.

How the Code Works

The C++, Python, and Java implementations all follow the same algorithm. They enumerate every increasing triple in the allowed range, discard non-primitive triples with a gcd test, and then evaluate the geometry for each surviving triple. For one triple, the implementation constructs the center triangle, derives the three tangency points, and applies the circumcenter formula to obtain $D$.

Next it applies Descartes' theorem to find the inner radius $r$, rewrites the three tangency conditions as two linear equations, and solves them to obtain $E$. The contribution of that triple is the Euclidean distance between $D$ and $E$. After all primitive triples have been processed, the accumulated sum is divided by the number of valid triples. One implementation also verifies several sample cases and checks that the loop count is $135739$.

Complexity Analysis

The search space consists of all triples with $1 \le a \lt b \lt c \le 100$, so the total enumeration cost is $O(100^3)$. For each triple, the program performs only a fixed number of arithmetic operations, square roots, and gcd evaluations, which is constant work. The memory usage is $O(1)$, since only a small set of coordinates and accumulators is kept at any time.

Footnotes and References

Problem page: https://projecteuler.net/problem=727
Descartes' circle theorem: Wikipedia - Descartes' theorem
Circumcenter and circumcircle formulas: Wikipedia - Circumscribed circle
Tangent circles: Wikipedia - Tangent circles
Euclidean distance: Wikipedia - Euclidean distance

Problem 727 source code

C++

#include <cassert>
#include <cstdint>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <algorithm>
#include <functional>

namespace {

constexpr double kEps = 1e-12;

struct Point {
    double x;
    double y;
};

double distance(const Point a, const Point b) {
    return std::hypot(a.x - b.x, a.y - b.y);
}

Point circumcenter(const Point p1, const Point p2, const Point p3) {
    const double x1 = p1.x;
    const double y1 = p1.y;
    const double x2 = p2.x;
    const double y2 = p2.y;
    const double x3 = p3.x;
    const double y3 = p3.y;

    const double den = 2.0 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2));

    const double n1 = x1 * x1 + y1 * y1;
    const double n2 = x2 * x2 + y2 * y2;
    const double n3 = x3 * x3 + y3 * y3;

    const double ux = (n1 * (y2 - y3) + n2 * (y3 - y1) + n3 * (y1 - y2)) / den;
    const double uy = (n1 * (x3 - x2) + n2 * (x1 - x3) + n3 * (x2 - x1)) / den;

    return {ux, uy};
}

double d_value(const int ra, const int rb, const int rc) {
    const double a = static_cast<double>(ra);
    const double b = static_cast<double>(rb);
    const double c = static_cast<double>(rc);

    const Point A{0.0, 0.0};
    const Point B{a + b, 0.0};

    const double cx = ((a + c) * (a + c) + (a + b) * (a + b) - (b + c) * (b + c)) / (2.0 * (a + b));
    const double cy = std::sqrt((a + c) * (a + c) - cx * cx);
    const Point C{cx, cy};

    const Point Pab{a, 0.0};
    const Point Pac{a / (a + c) * cx, a / (a + c) * cy};
    const Point Pbc{B.x + b / (b + c) * (cx - B.x), b / (b + c) * cy};

    const Point D = circumcenter(Pab, Pac, Pbc);

    const double k1 = 1.0 / a;
    const double k2 = 1.0 / b;
    const double k3 = 1.0 / c;
    const double k4 = k1 + k2 + k3 + 2.0 * std::sqrt(k1 * k2 + k2 * k3 + k3 * k1);
    const double r = 1.0 / k4;

    const double rhs_ab = (B.x * B.x + B.y * B.y) + (a + r) * (a + r) - (b + r) * (b + r);
    const double ex = rhs_ab / (2.0 * B.x);

    const double rhs_ac = (C.x * C.x + C.y * C.y) + (a + r) * (a + r) - (c + r) * (c + r);
    const double ey = (rhs_ac - 2.0 * C.x * ex) / (2.0 * C.y);
    const Point E{ex, ey};

    assert(std::fabs(distance(E, A) - (a + r)) < 1e-9);
    assert(std::fabs(distance(E, B) - (b + r)) < 1e-9);
    assert(std::fabs(distance(E, C) - (c + r)) < 1e-9);

    return distance(D, E);
}

double expected_d() {
    double sum = 0.0;
    std::int64_t count = 0;

    for (int a = 1; a <= 100; ++a) {
        for (int b = a + 1; b <= 100; ++b) {
            for (int c = b + 1; c <= 100; ++c) {
                if (std::gcd(a, std::gcd(b, c)) != 1) {
                    continue;
                }
                sum += d_value(a, b, c);
                ++count;
            }
        }
    }

    assert(count == 135739);
    return sum / static_cast<double>(count);
}

}  // namespace

int main() {
    assert(std::fabs(d_value(1, 2, 3) - 0.15676309893321677) < 1e-12);
    assert(std::fabs(d_value(1, 2, 4) - 0.19528865932762782) < 1e-12);
    assert(std::fabs(d_value(2, 3, 5) - 0.2060868349032223) < 1e-12);

    std::cout << std::fixed << std::setprecision(8) << expected_d() << '\n';
    return 0;
}

Python

import math

class Point:
    def __init__(self, x, y):
        self.x = x
        self.y = y

def distance(a, b):
    return math.hypot(a.x - b.x, a.y - b.y)

def circumcenter(p1, p2, p3):
    x1, y1 = p1.x, p1.y
    x2, y2 = p2.x, p2.y
    x3, y3 = p3.x, p3.y
    
    den = 2.0 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2))
    
    n1 = x1 * x1 + y1 * y1
    n2 = x2 * x2 + y2 * y2
    n3 = x3 * x3 + y3 * y3
    
    ux = (n1 * (y2 - y3) + n2 * (y3 - y1) + n3 * (y1 - y2)) / den
    uy = (n1 * (x3 - x2) + n2 * (x1 - x3) + n3 * (x2 - x1)) / den
    
    return Point(ux, uy)

def d_value(ra, rb, rc):
    a = float(ra)
    b = float(rb)
    c = float(rc)
    
    A = Point(0.0, 0.0)
    B = Point(a + b, 0.0)
    
    cx = ((a + c) * (a + c) + (a + b) * (a + b) - (b + c) * (b + c)) / (2.0 * (a + b))
    cy = math.sqrt(max(0.0, (a + c) * (a + c) - cx * cx))
    C = Point(cx, cy)
    
    Pab = Point(a, 0.0)
    Pac = Point(a / (a + c) * cx, a / (a + c) * cy)
    Pbc = Point(B.x + b / (b + c) * (cx - B.x), b / (b + c) * cy)
    
    D = circumcenter(Pab, Pac, Pbc)
    
    k1 = 1.0 / a
    k2 = 1.0 / b
    k3 = 1.0 / c
    k4 = k1 + k2 + k3 + 2.0 * math.sqrt(k1 * k2 + k2 * k3 + k3 * k1)
    r = 1.0 / k4
    
    rhs_ab = (B.x * B.x + B.y * B.y) + (a + r) * (a + r) - (b + r) * (b + r)
    ex = rhs_ab / (2.0 * B.x)
    
    rhs_ac = (C.x * C.x + C.y * C.y) + (a + r) * (a + r) - (c + r) * (c + r)
    ey = (rhs_ac - 2.0 * C.x * ex) / (2.0 * C.y)
    E = Point(ex, ey)
    
    return distance(D, E)

def solve():
    total_sum = 0.0
    count = 0
    
    for a in range(1, 101):
        for b in range(a + 1, 101):
            for c in range(b + 1, 101):
                if math.gcd(a, math.gcd(b, c)) != 1:
                    continue
                total_sum += d_value(a, b, c)
                count += 1
                
    ans = total_sum / count
    return f"{ans:.8f}"

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

Java

public class Euler727 {
    static class Point {
        double x, y;

        Point(double x, double y) {
            this.x = x;
            this.y = y;
        }
    }

    static double distance(Point a, Point b) {
        return Math.hypot(a.x - b.x, a.y - b.y);
    }

    static Point circumcenter(Point p1, Point p2, Point p3) {
        double x1 = p1.x, y1 = p1.y;
        double x2 = p2.x, y2 = p2.y;
        double x3 = p3.x, y3 = p3.y;

        double den = 2.0 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2));

        double n1 = x1 * x1 + y1 * y1;
        double n2 = x2 * x2 + y2 * y2;
        double n3 = x3 * x3 + y3 * y3;

        double ux = (n1 * (y2 - y3) + n2 * (y3 - y1) + n3 * (y1 - y2)) / den;
        double uy = (n1 * (x3 - x2) + n2 * (x1 - x3) + n3 * (x2 - x1)) / den;

        return new Point(ux, uy);
    }

    static double dValue(int ra, int rb, int rc) {
        double a = ra, b = rb, c = rc;

        Point A = new Point(0.0, 0.0);
        Point B = new Point(a + b, 0.0);

        double cx = ((a + c) * (a + c) + (a + b) * (a + b) - (b + c) * (b + c)) / (2.0 * (a + b));
        double cy = Math.sqrt(Math.max(0.0, (a + c) * (a + c) - cx * cx));
        Point C = new Point(cx, cy);

        Point Pab = new Point(a, 0.0);
        Point Pac = new Point(a / (a + c) * cx, a / (a + c) * cy);
        Point Pbc = new Point(B.x + b / (b + c) * (cx - B.x), b / (b + c) * cy);

        Point D = circumcenter(Pab, Pac, Pbc);

        double k1 = 1.0 / a;
        double k2 = 1.0 / b;
        double k3 = 1.0 / c;
        double k4 = k1 + k2 + k3 + 2.0 * Math.sqrt(k1 * k2 + k2 * k3 + k3 * k1);
        double r = 1.0 / k4;

        double rhs_ab = (B.x * B.x + B.y * B.y) + (a + r) * (a + r) - (b + r) * (b + r);
        double ex = rhs_ab / (2.0 * B.x);

        double rhs_ac = (C.x * C.x + C.y * C.y) + (a + r) * (a + r) - (c + r) * (c + r);
        double ey = (rhs_ac - 2.0 * C.x * ex) / (2.0 * C.y);
        Point E = new Point(ex, ey);

        return distance(D, E);
    }

    static long gcd(long a, long b) {
        while (b != 0) {
            long temp = b;
            b = a % b;
            a = temp;
        }
        return a;
    }

    public static String solve() {
        double sum = 0.0;
        long count = 0;

        for (int a = 1; a <= 100; ++a) {
            for (int b = a + 1; b <= 100; ++b) {
                for (int c = b + 1; c <= 100; ++c) {
                    if (gcd(a, gcd(b, c)) != 1) {
                        continue;
                    }
                    sum += dValue(a, b, c);
                    ++count;
                }
            }
        }

        double ans = sum / count;
        return String.format(java.util.Locale.US, "%.8f", ans);
    }

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