Solution

Answer: 315306518862563689

Let

$r_A=a$,
$r_B=b$,
$r_C=c$,

with $0<a\le b\le n$, and all radii integers.

We must compute

$$S(n)=\sum (a+b+c)$$

over all valid tangent configurations.

1. Mathematical analysis

We place the tangent line $L$ as the $x$-axis.

Then:

circle $A$ has center $(0,a)$,
circle $B$ has center $(2\sqrt{ab},,b)$,

because two circles tangent to the same line and to each other satisfy:

$$(\Delta x)^2=(a+b)^2-(a-b)^2=4ab.$$

Now place circle $C$ with center $(x,c)$.

Since $C$ is tangent to $A$,

$$x^2=(a+c)^2-(a-c)^2=4ac,$$

hence

$$x=2\sqrt{ac}.$$

Likewise, tangency with $B$ gives

$$(2\sqrt{ab}-x)^2=4bc.$$

Substitute $x=2\sqrt{ac}$:

$$(2\sqrt{ab}-2\sqrt{ac})^2=4bc.$$

Divide by $4$:

$$(\sqrt{ab}-\sqrt{ac})^2=bc.$$

Taking square roots in the geometric configuration yields

$$\sqrt{ab}=\sqrt{ac}+\sqrt{bc}.$$

Factor out $\sqrt c$:

$$\sqrt{ab}=\sqrt c(\sqrt a+\sqrt b).$$

Therefore

$$\boxed{\sqrt c=\frac{\sqrt{ab}}{\sqrt a+\sqrt b}}$$

and

$$\boxed{c=\frac{ab}{(\sqrt a+\sqrt b)^2}}.$$

Integer parametrization

For $c$ to be integral, write

$$a=dx^2,\qquad b=dy^2,$$

where $x,y$ are coprime.

Then

$$c=\frac{d^2x^2y^2}{d(x+y)^2} =\frac{dx^2y^2}{(x+y)^2}.$$

Since

$$\gcd(x+y,xy)=1$$

for coprime $x,y$, we must have

$$(x+y)^2\mid d.$$

So let

$$d=k(x+y)^2.$$

Then every solution is

$$\boxed{ \begin{aligned} a&=k,x^2(x+y)^2,\ b&=k,y^2(x+y)^2,\ c&=k,x^2y^2, \end{aligned}}$$

with

$$\gcd(x,y)=1,\qquad x\le y,\qquad k\ge1.$$

Contribution to $S(n)$

Since $b\le n$,

$$k,y^2(x+y)^2\le n.$$

Define

$$M=\left\lfloor \frac{n}{y^2(x+y)^2}\right\rfloor.$$

For fixed coprime $x,y$, the valid values are $k=1,\dots,M$.

Now

$$a+b+c = k\Big((x+y)^2(x^2+y^2)+x^2y^2\Big).$$

Hence the total contribution is

$$\Big((x+y)^2(x^2+y^2)+x^2y^2\Big) \sum_{k=1}^{M}k.$$

Using

$$\sum_{k=1}^{M}k=\frac{M(M+1)}2,$$

we obtain

$$\boxed{ S(n)= \sum_{\substack{x\le y\ \gcd(x,y)=1}} \left((x+y)^2(x^2+y^2)+x^2y^2\right) \frac{M(M+1)}2 }$$

where

$$M=\left\lfloor \frac{n}{y^2(x+y)^2}\right\rfloor.$$

This is efficient because

$$y^2(x+y)^2\le n$$

implies only about $O(\sqrt n)$ relevant pairs.

2. Python implementation

from math import gcd, isqrt

def S(n):
    total = 0

    # y^2 (x+y)^2 <= n
    # so y is at most sqrt(n)
    limit = isqrt(n)

    for y in range(1, limit + 1):
        for x in range(1, y + 1):

            # primitive pairs only
            if gcd(x, y) != 1:
                continue

            denom = y * y * (x + y) * (x + y)

            if denom > n:
                break

            M = n // denom

            base = (x + y) * (x + y) * (x * x + y * y)
            base += x * x * y * y

            total += base * M * (M + 1) // 2

    return total

print(S(5))       # 9
print(S(100))     # 3072
print(S(10**9))

3. Code walkthrough

Imports

from math import gcd, isqrt

gcd checks coprimality of $x,y$.
isqrt computes integer square roots efficiently.

Main function

def S(n):
    total = 0

Accumulates the answer.

Bound for iteration

limit = isqrt(n)

Because

$$y^2(x+y)^2 \le n,$$

we certainly have $y \le \sqrt n$.

Enumerating primitive pairs

for y in range(1, limit + 1):
    for x in range(1, y + 1):

We enforce $x\le y$, corresponding to $a\le b$.

Coprime condition

if gcd(x, y) != 1:
    continue

Required by the parametrization.

Denominator constraint

denom = y * y * (x + y) * (x + y)

This equals

$$y^2(x+y)^2.$$

If it exceeds $n$, no scaling factor $k\ge1$ works.

Maximum scaling factor

M = n // denom

This is

$$\left\lfloor \frac{n}{y^2(x+y)^2}\right\rfloor.$$

Base contribution

base = (x + y) * (x + y) * (x * x + y * y)
base += x * x * y * y

This equals

$$(x+y)^2(x^2+y^2)+x^2y^2.$$

Sum over all $k$

total += base * M * (M + 1) // 2

Using

$$1+2+\cdots+M=\frac{M(M+1)}2.$$

4. Verification

The program reproduces the given values:

$$S(5)=9$$

and

$$S(100)=3072.$$

So the derivation and implementation are consistent.

Running the computation for $n=10^9$ gives:

$$315306518862563689.$$

The magnitude is reasonable since the summation grows roughly quadratically with scaling contributions.

Answer: 315306518862563689

Project Euler Problem 510