Project Euler Problem 727

Project Euler Problem 727

Solution

Answer: 3.64039141

Let the three mutually tangent circles have centres $A,B,C$ and radii $a,b,c$ respectively.

We place the circles in coordinates as follows:

• $A=(0,0)$
• $B=(a+b,0)$
• $C=(x,y)$

Since the circles are externally tangent,

$$|AC|=a+c,\qquad |BC|=b+c.$$

So

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

$$y=\sqrt{(a+c)^2-x^2}.$$

The tangency points are:

• between $A,B$:

$$P_{ab}=(a,0)$$

• between $A,C$:

$$P_{ac}=\frac{a}{a+c}C$$

• between $B,C$:

$$P_{bc}=B+\frac{b}{b+c}(C-B).$$

The circumcenter $D$ is the circumcenter of triangle

$(P_{ab},P_{ac},P_{bc})$.

The inner Soddy circle radius is given by Descartes’ theorem:

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

where $k_i=1/r_i$, and therefore

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

Its center $E$ satisfies

$$|EA|=a+r,\quad |EB|=b+r.$$

Finally,

$$d=|DE|.$$

We average $d$ over all coprime triples

$$1\le a<b<c\le100.$$

Python code used:

from math import sqrt, gcd

def dist(p, q):
    return ((p[0]-q[0])**2 + (p[1]-q[1])**2) ** 0.5

def circumcenter(p1, p2, p3):
    ax, ay = p1
    bx, by = p2
    cx, cy = p3

    d = 2 * (ax*(by-cy) + bx*(cy-ay) + cx*(ay-by))

    ux = (
        (ax*ax + ay*ay)*(by-cy)
        + (bx*bx + by*by)*(cy-ay)
        + (cx*cx + cy*cy)*(ay-by)
    ) / d

    uy = (
        (ax*ax + ay*ay)*(cx-bx)
        + (bx*bx + by*by)*(ax-cx)
        + (cx*cx + cy*cy)*(bx-ax)
    ) / d

    return (ux, uy)

def compute(a, b, c):
    A = (0.0, 0.0)
    B = (a + b, 0.0)

    x = ((a+c)**2 + (a+b)**2 - (b+c)**2) / (2*(a+b))
    y = sqrt((a+c)**2 - x*x)

    C = (x, y)

    Pab = (a, 0.0)

    Pac = (
        a*x/(a+c),
        a*y/(a+c)
    )

    Pbc = (
        B[0] + b*(x-B[0])/(b+c),
        b*y/(b+c)
    )

    D = circumcenter(Pab, Pac, Pbc)

    k1, k2, k3 = 1/a, 1/b, 1/c

    k4 = (
        k1 + k2 + k3
        + 2*sqrt(k1*k2 + k2*k3 + k3*k1)
    )

    r = 1 / k4

    ex = (
        (a+r)**2 - (b+r)**2 + (a+b)**2
    ) / (2*(a+b))

    ey = sqrt((a+r)**2 - ex*ex)

    E = (ex, ey)

    return dist(D, E)

total = 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 gcd(gcd(a, b), c) == 1:
                total += compute(a, b, c)
                count += 1

answer = total / count
print(answer)

This produces

$$3.640391406429808$$

which rounds to eight decimal places as:

Answer: 3.64039141