Problem 510: Tangent Circles
View on Project EulerProject Euler Problem 510 Solution
EulerSolve provides an optimized solution for Project Euler Problem 510, Tangent Circles, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Let \(S(N)\) denote the sum of \(a+b+c\) over all integer triples \((a,b,c)\) with \(1\le a\le b\le N\) such that $$\sqrt{ab}\in \mathbb{Z},\qquad c=\frac{ab}{a+b+2\sqrt{ab}}\in \mathbb{Z}_{>0}.$$ The tangent-circle geometry is already condensed into this arithmetic condition. The difficulty is that \(a\) and \(b\) are tied together by both a square-root constraint and a divisibility constraint, so scanning all pairs up to \(N=10^9\) is far too slow. The key is to show that every valid triple comes from one coprime parameter pair and one scale factor. Mathematical Approach The solution derives a complete parameterization and then sums one whole family at a time. Step 1: Separate the Square Part of \(ab\) Write $$a=g u,\qquad b=g v,\qquad \gcd(u,v)=1,$$ where \(g=\gcd(a,b)\). Because \(\sqrt{ab}\) is an integer, the product \(ab=g^2uv\) is a perfect square, so \(uv\) is also a perfect square....
Detailed mathematical approach
Problem Summary
Let \(S(N)\) denote the sum of \(a+b+c\) over all integer triples \((a,b,c)\) with \(1\le a\le b\le N\) such that
$$\sqrt{ab}\in \mathbb{Z},\qquad c=\frac{ab}{a+b+2\sqrt{ab}}\in \mathbb{Z}_{>0}.$$
The tangent-circle geometry is already condensed into this arithmetic condition. The difficulty is that \(a\) and \(b\) are tied together by both a square-root constraint and a divisibility constraint, so scanning all pairs up to \(N=10^9\) is far too slow. The key is to show that every valid triple comes from one coprime parameter pair and one scale factor.
Mathematical Approach
The solution derives a complete parameterization and then sums one whole family at a time.
Step 1: Separate the Square Part of \(ab\)
Write
$$a=g u,\qquad b=g v,\qquad \gcd(u,v)=1,$$
where \(g=\gcd(a,b)\). Because \(\sqrt{ab}\) is an integer, the product \(ab=g^2uv\) is a perfect square, so \(uv\) is also a perfect square. Since \(u\) and \(v\) are coprime, each of them must already be a square:
$$u=m^2,\qquad v=n^2,\qquad \gcd(m,n)=1.$$
Therefore every valid pair can be written as
$$a=g m^2,\qquad b=g n^2,\qquad \sqrt{ab}=gmn.$$
Step 2: Force the Denominator to Divide Exactly
Substituting the previous form into the definition of \(c\) gives
$$c=\frac{g^2m^2n^2}{g(m^2+n^2+2mn)}=\frac{g m^2n^2}{(m+n)^2}.$$
Now
$$\gcd(m+n,m)=\gcd(m+n,n)=1,$$
so \((m+n)\) is coprime to both \(m\) and \(n\), hence
$$\gcd\left((m+n)^2,m^2n^2\right)=1.$$
That means the entire denominator \((m+n)^2\) must come from the common factor \(g\). So there is a positive integer \(t\) with
$$g=t(m+n)^2.$$
We obtain the exact shape of every valid triple:
$$a=t(m+n)^2m^2,\qquad b=t(m+n)^2n^2,\qquad c=t m^2n^2.$$
Step 3: Show the Parameterization Is Complete and Unique
The converse is immediate. If \(\gcd(m,n)=1\) and \(t\ge 1\), then
$$\sqrt{ab}=t(m+n)^2mn\in \mathbb{Z},$$
and
$$a+b+2\sqrt{ab}=t(m+n)^2(m^2+n^2+2mn)=t(m+n)^4.$$
Therefore
$$\frac{ab}{a+b+2\sqrt{ab}}=\frac{t^2(m+n)^4m^2n^2}{t(m+n)^4}=t m^2n^2=c,$$
so every triple of this form is valid. The coprime pair \((m,n)\) is unique up to swapping the two entries, so enforcing
$$m\le n$$
removes double counting. For one coprime pair, the primitive triple is
$$a_0=(m+n)^2m^2,\qquad b_0=(m+n)^2n^2,\qquad c_0=m^2n^2.$$
Step 4: Reduce the Bound to One Variable
When \(m\le n\), the three components satisfy
$$c_0\le a_0\le b_0.$$
So after scaling by \(t\), the condition that all components lie under \(N\) reduces to the single inequality
$$t b_0\le N.$$
Hence the admissible scales are exactly
$$1\le t\le t_{\max}=\left\lfloor\frac{N}{b_0}\right\rfloor.$$
The entire family generated by one primitive triple contributes
$$\sum_{t=1}^{t_{\max}} t(a_0+b_0+c_0)=(a_0+b_0+c_0)\frac{t_{\max}(t_{\max}+1)}{2}.$$
Step 5: Enumerate Only the Feasible Coprime Pairs
For fixed \(n\), the quantity \(b_0=(m+n)^2n^2\) increases with \(m\), so once \(b_0>N\) the inner search can stop immediately. The smallest possible \(b_0\) for that \(n\) occurs at \(m=1\):
$$b_{0,\min}=(n+1)^2n^2.$$
If this already exceeds \(N\), then no larger \(n\) can work. Because \(b_{0,\min}\) is on the order of \(n^4\), only values
$$n=O(N^{1/4})$$
must be inspected.
Worked Example: \(N=100\)
The feasible coprime pairs with \(m\le n\) are very few.
For \((m,n)=(1,1)\), the primitive triple is
$$(a_0,b_0,c_0)=(4,4,1),\qquad t_{\max}=\left\lfloor\frac{100}{4}\right\rfloor=25.$$
Its contribution is
$$(4+4+1)\frac{25\cdot 26}{2}=9\cdot 325=2925.$$
For \((m,n)=(1,2)\), the primitive triple is
$$(a_0,b_0,c_0)=(9,36,4),\qquad t_{\max}=\left\lfloor\frac{100}{36}\right\rfloor=2,$$
so the contribution is
$$(9+36+4)\frac{2\cdot 3}{2}=49\cdot 3=147.$$
The pair \((2,2)\) is excluded because it is not coprime, and for \(n=3\) even the smallest possible largest term is
$$(3+1)^2\cdot 3^2=144>100.$$
Therefore
$$S(100)=2925+147=3072,$$
matching the checkpoint used by the implementation.
How the Code Works
The C++, Python, and Java implementations all use the parameterization above instead of testing arbitrary pairs \((a,b)\). They iterate the larger coprime parameter upward, stop as soon as \((n+1)^2n^2>N\), and for each admissible value scan the smaller parameter up to \(n\).
Non-coprime pairs are skipped immediately. For each remaining pair, the implementation builds the primitive triple \((a_0,b_0,c_0)\), breaks the inner loop when \(b_0>N\), computes
$$t_{\max}=\left\lfloor\frac{N}{b_0}\right\rfloor,$$
and adds
$$(a_0+b_0+c_0)\frac{t_{\max}(t_{\max}+1)}{2}.$$
So the code never searches invalid triples and never performs an expensive global brute-force scan. It sums complete scaling families directly.
Complexity Analysis
Let \(R=\lfloor N^{1/4}\rfloor\). Because the search only considers \(1\le m\le n\le R\), the number of candidate pairs is at most \(R(R+1)/2=O(\sqrt{N})\). Each pair needs one coprimality test and a constant amount of integer arithmetic, so the overall running time is \(O(\sqrt{N})\) in arithmetic operations, with \(O(1)\) extra memory.
Footnotes and References
- Problem page: https://projecteuler.net/problem=510
- Greatest common divisor: Wikipedia — Greatest common divisor
- Coprime integers: Wikipedia — Coprime integers
- Square number: Wikipedia — Square number
- Triangular number: Wikipedia — Triangular number
Problem 510 source code
C++
#include <cmath>
#include <cstdint>
#include <iostream>
#include <numeric>
namespace {
using u64 = std::uint64_t;
using u128 = __uint128_t;
u64 isqrt_u128(const u128 x) {
long double approx = std::sqrt(static_cast<long double>(x));
u64 r = static_cast<u64>(approx);
while (static_cast<u128>(r + 1ULL) * static_cast<u128>(r + 1ULL) <= x) {
++r;
}
while (static_cast<u128>(r) * static_cast<u128>(r) > x) {
--r;
}
return r;
}
u64 sum_formula(const u64 n_max) {
u128 total = 0;
for (u64 n = 1ULL;; ++n) {
const u128 min_b = static_cast<u128>(n + 1ULL) * static_cast<u128>(n + 1ULL) *
static_cast<u128>(n) * static_cast<u128>(n);
if (min_b > n_max) {
break;
}
for (u64 m = 1ULL; m <= n; ++m) {
if (std::gcd(m, n) != 1ULL) {
continue;
}
const u128 sum = static_cast<u128>(m + n);
const u128 b0 = sum * sum * static_cast<u128>(n) * static_cast<u128>(n);
if (b0 > n_max) {
break;
}
const u128 a0 = sum * sum * static_cast<u128>(m) * static_cast<u128>(m);
const u128 c0 = static_cast<u128>(m) * static_cast<u128>(m) * static_cast<u128>(n) *
static_cast<u128>(n);
const u128 t_max = static_cast<u128>(n_max) / b0;
const u128 t_sum = t_max * (t_max + 1ULL) / 2ULL;
total += (a0 + b0 + c0) * t_sum;
}
}
return static_cast<u64>(total);
}
u64 sum_bruteforce(const u64 n_max) {
u128 total = 0;
for (u64 a = 1ULL; a <= n_max; ++a) {
for (u64 b = a; b <= n_max; ++b) {
const u128 prod = static_cast<u128>(a) * static_cast<u128>(b);
const u64 s = isqrt_u128(prod);
if (static_cast<u128>(s) * static_cast<u128>(s) != prod) {
continue;
}
const u128 denom = static_cast<u128>(a + b) + 2ULL * static_cast<u128>(s);
if (prod % denom != 0ULL) {
continue;
}
const u128 c = prod / denom;
total += static_cast<u128>(a + b) + c;
}
}
return static_cast<u64>(total);
}
bool run_checkpoints() {
if (sum_formula(5ULL) != 9ULL) {
std::cerr << "Checkpoint failed: S(5)\n";
return false;
}
if (sum_formula(100ULL) != 3'072ULL) {
std::cerr << "Checkpoint failed: S(100)\n";
return false;
}
if (sum_formula(300ULL) != sum_bruteforce(300ULL)) {
std::cerr << "Checkpoint failed: formula/bruteforce mismatch at 300\n";
return false;
}
return true;
}
} // namespace
int main() {
if (!run_checkpoints()) {
return 1;
}
constexpr u64 n = 1'000'000'000ULL;
std::cout << sum_formula(n) << '\n';
return 0;
}
Python
import math
def solve():
N_MAX = 1_000_000_000
total = 0
n = 1
while True:
min_b = (n + 1)**2 * n**2
if min_b > N_MAX:
break
for m in range(1, n + 1):
if math.gcd(m, n) != 1:
continue
s = m + n
b0 = s * s * n * n
if b0 > N_MAX:
break
a0 = s * s * m * m
c0 = m * m * n * n
t_max = N_MAX // b0
t_sum = t_max * (t_max + 1) // 2
total += (a0 + b0 + c0) * t_sum
n += 1
return str(total)
if __name__ == '__main__':
print(solve())
Java
public class Euler510 {
private static long gcd(long a, long b) {
while (b != 0) {
long t = b;
b = a % b;
a = t;
}
return a;
}
private static long sumFormula(long nMax) {
long total = 0;
for (long n = 1;; ++n) {
long minB = (n + 1) * (n + 1) * n * n;
if (minB > nMax)
break;
for (long m = 1; m <= n; ++m) {
if (gcd(m, n) != 1)
continue;
long s = m + n;
long b0 = s * s * n * n;
if (b0 > nMax)
break;
long a0 = s * s * m * m;
long c0 = m * m * n * n;
long tMax = nMax / b0;
long tSum = tMax * (tMax + 1) / 2;
total += (a0 + b0 + c0) * tSum;
}
}
return total;
}
public static void main(String[] args) {
long nMax = 1000000000L;
System.out.println(sumFormula(nMax));
}
}