Problem 373: Circumscribed Circles

Project Euler Problem 373 Solution

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

Problem Summary For each integer circumradius $r$, let $T(r)$ be the number of triangles with integer side lengths and circumradius exactly $r$. The program computes $$S(N)=\sum_{r=1}^{N} r\,T(r),\qquad N=10^7.$$ A brute-force search over all triples $(a,b,c)$ up to $2r$ would be far too slow, so the implementation first characterizes which side lengths can occur for a fixed radius and then checks only those candidates exactly. Mathematical Approach Step 1: Convert the circumradius condition into an exact integer identity For a triangle with sides $a,b,c$ and area $\Delta$, Heron's formula gives $$16\Delta^2=(a+b+c)(a+b-c)(a-b+c)(-a+b+c).$$ The circumradius also satisfies $$R=\frac{abc}{4\Delta}.$$ Squaring and eliminating $\Delta$ yields $$a^2b^2c^2=R^2(a+b+c)(a+b-c)(a-b+c)(-a+b+c).$$ This is exactly the predicate used in the solver. It avoids floating-point roundoff entirely: once $a$, $b$, $c$, and $r$ are integers, the question “is the circumradius equal to $r$?” becomes a pure integer comparison. Step 2: Describe every possible side length for a fixed radius If side $a$ is opposite angle $A$, then the extended law of sines gives $$a=2r\sin A.$$ Because $a$ and $r$ are integers, $\sin A=a/(2r)$ is rational....

Detailed mathematical approach

Problem Summary

Every triangle has exactly one circumscribed circle through its three vertices; its radius is the circumradius. Among all triangles whose three side lengths are integers we keep those whose circumradius is an integer as well, and we let $S(n)$ be the sum of the radii of all of them whose radius does not exceed $n$. Triangles are counted up to congruence, so a triangle is just an unordered triple of side lengths. The statement gives the checkpoints $S(100)=4950$ and $S(1200)=1653605$, and asks for $S(10^7)$. ^[1]

Writing $T(r)$ for the number of integer-sided triangles whose circumradius is exactly $r$, the quantity we need is

$$S(n)=\sum_{r=1}^{n} r\,T(r).$$

The whole solution in one sentence: for a fixed radius $r$ an integer side can only be the diameter $2r$ or $2ru/d$, where $u$ and $d$ are a leg and the hypotenuse of a primitive Pythagorean triple with $d$ dividing $r$; that list of candidates is very short, so we build it and test each of its triples with one exact identity between integers.

Mathematical Approach

1. An exact test for the circumradius

Floating-point trigonometry cannot decide whether a radius is exactly an integer, so the first task is to turn the question into a comparison of integers. Heron's formula writes the area $\Delta$ of a triangle with sides $a,b,c$ as ^[3]

$$\begin{aligned}16\Delta^2={}&(a+b+c)(a+b-c)\\&\cdot(a-b+c)(-a+b+c),\end{aligned}$$

and the circumradius of the same triangle satisfies ^[2]

$$R=\frac{abc}{4\Delta}.$$

Squaring the second identity and substituting the first eliminates the area and leaves a statement about integers only:

$$\begin{aligned}a^2b^2c^2=R^2{}&(a+b+c)(a+b-c)\\&\cdot(a-b+c)(-a+b+c).\end{aligned}$$

For the right triangle $a,b,c=6,8,10$ and $R=5$ the left side is $36\cdot64\cdot100=230400$ and the right side is $25\cdot24\cdot4\cdot8\cdot12=230400$, so the test accepts it. The four bracketed factors are positive exactly when $a,b,c$ obey the strict triangle inequality, ^[3] and for such sides the area is a fixed positive number, so the identity holds for exactly one positive $R$. So once the triangle inequality is checked, the identity with $R=r$ holds if and only if the circumradius is $r$. This test is exact but far too slow for all triples of sides up to $2\cdot10^7$; the rest of the mathematics shrinks the set it must see.

2. Every side is a chord, so its sine is rational

Place the triangle in its circumscribed circle of radius $r$. Each side is a chord, and the extended law of sines relates it to the angle opposite and to the radius: ^[4]

$$\begin{aligned}a&=2r\sin A,\quad b=2r\sin B,\\ c&=2r\sin C.\end{aligned}$$

If $a$ and $r$ are integers then $\sin A=a/(2r)$ is a rational number. Since $0\lt A\lt\pi$ we have $0\lt\sin A\le1$, and $\sin A=1$ means $A=\pi/2$, in which case the side is the full diameter $a=2r$.

3. The cosine is rational too

This is the step that cuts the candidates down from all integers up to $2r$ to a mere handful. The law of cosines gives ^[5]

$$\cos A=\frac{b^2+c^2-a^2}{2bc},$$

whose right-hand side is a quotient of integers. So for every angle both the sine and the cosine are rational: $(\cos A,\sin A)$ is a rational point on the unit circle.

4. A rational point on the circle is a primitive Pythagorean triple

Write $\sin A=p/q$ and $\cos A=s/q'$ in lowest terms, with $q,q'\gt0$ and $s$ possibly negative or zero. We claim $q=q'$.

Indeed $p^2/q^2+s^2/q'^2=1$ gives $p^2q'^2+s^2q^2=q^2q'^2$. Reading it modulo $q^2$ shows $q^2\mid p^2q'^2$, and $\gcd(p,q)=1$ forces $q^2\mid q'^2$; the symmetric argument gives $q'^2\mid q^2$, so $q^2=q'^2$ and $q=q'$. Cancelling $q^2$ then leaves

$$p^2+s^2=q^2.$$

If $s=0$ then $p=q=1$, the diameter case. Otherwise $p$ and $|s|$ are the legs and $q$ the hypotenuse of a Pythagorean triple, that is of a right triangle with integer sides; the triple is primitive, its entries sharing no divisor above $1$, because $\gcd(p,q)=1$. So apart from the diameter, such a sine is always a leg divided by the hypotenuse of a primitive triple.

5. The legs belonging to one hypotenuse

Euclid's parametrisation is complete: for integers $m\gt n\ge1$ that are coprime, so $\gcd(m,n)=1$, and of opposite parity, one even and one odd, ^[6]

$$\begin{aligned}u_1&=m^2-n^2,\quad u_2=2mn,\\ d&=m^2+n^2\end{aligned}$$

runs through every primitive triple, with the two legs $u_1,u_2$ in either role. For $(m,n)=(2,1)$ this is $3,4,5$ and for $(m,n)=(4,3)$ it is $7,24,25$. A given $d$ may carry several such pairs: $65=8^2+1^2=7^2+4^2$ yields $63,16,65$ and $33,56,65$. The number of primitive triples with a fixed hypotenuse $d$ is $2^{k-1}$, where $k$ is the number of distinct prime factors of $d$, ^[6] so $d$ carries exactly $2^{k}$ legs.

6. The hypotenuse is odd, so it must divide the radius

A primitive hypotenuse is always odd. ^[6] The legs cannot both be even, by primitivity, nor both odd, since then $q^2=p^2+s^2\equiv2\pmod 4$, which no square is; so one is even and one odd, and $q^2$ is odd.

Now suppose a side satisfies $a=2r\,u/d$ with $\gcd(u,d)=1$. For $a$ to be an integer we need $d\mid 2ru$, hence $d\mid 2r$, and because $d$ is odd this is the same as

$$d\mid r,\qquad a=2\frac{r}{d}\,u.$$

That is the decisive restriction: only the divisors of $r$ need inspecting, not all hypotenuses up to $r$.

7. The candidate list for one radius

Call $d$ a primitive hypotenuse if it is the hypotenuse of some primitive triple. These are exactly the integers $d\gt1$ all of whose prime factors are congruent to $1$ modulo $4$: ^[6] ^[7] they begin $5,13,17,25,29,37,41,53,61,65$. Collecting Sections 2 to 6, every side of an integer-sided triangle with circumradius $r$ lies in

$$\mathcal S_r=\{2r\}\cup\Bigl\{2\tfrac{r}{d}u\ :\ d\mid r,\ u\in U(d)\Bigr\},$$

where $U(d)$ is the set of legs of the primitive triples with hypotenuse $d$, and is empty when $d$ is not a primitive hypotenuse.

These candidates are automatically distinct: each is $2r$ times a fraction $u/d$ already in lowest terms, a leg and the hypotenuse of a primitive triple being coprime, and the diameter is the case $u/d=1/1$, so two equal candidates share the same $d$ and $u$. Counting them with the $2^{k}$ legs of Section 5, and letting $p^{e}$ run over the prime powers of $r$ with $p\equiv1\pmod 4$,

$$|\mathcal S_r|=\prod_{p\equiv1\ (4)}(1+2e).$$

Each such prime offers the exponent $0$, giving the factor $1$, or one of $e$ positive exponents, each doubling the count. For example $|\mathcal S_{25}|=5$, $|\mathcal S_{1105}|=27$ because $1105=5\cdot13\cdot17$, and the largest value for $r\le10^7$ is $405$, reached at $r=5928325$, which is $5^2\cdot13\cdot17\cdot29\cdot37$.

8. A worked example: $r=25$

The divisors of $25$ are $1,5,25$; of these $5$ and $25$ are primitive hypotenuses. Divisor $5$ comes from $(m,n)=(2,1)$ with legs $3,4$, and the scale factor is $2\cdot(25/5)=10$, giving the sides $30$ and $40$. Divisor $25$ comes from $(m,n)=(4,3)$ with legs $7,24$ and scale factor $2$, giving $14$ and $48$. With the diameter $50$ the candidate list is

$$\mathcal S_{25}=\{14,30,40,48,50\}.$$

Repeated sides are allowed, so the five candidates give $\binom{5+2}{3}=35$ unordered triples: sending $i\le j\le k$ to $i,j+1,k+2$ pairs them with the three-element subsets of $1,\dots,7$. Six of them fail the triangle inequality, for instance $14+14\le30$; the remaining $29$ go to the exact test of Section 1, and exactly five pass:

$(14,48,50)$ and $(30,40,50)$, right-angled on the diameter, since $14^2+48^2=50^2$ and $30^2+40^2=50^2$;
$(30,30,48)$, isosceles and obtuse, and $(40,40,48)$, isosceles and acute;
$(14,30,40)$, obtuse, with $\cos$ of its largest angle equal to $-3/5$.

So $T(25)=5$, which contributes $125$ to $S(n)$ for every $n\ge25$. The exact test cannot be skipped: $(14,14,14)$ obeys the triangle inequality but has circumradius $14/\sqrt3$, so it is rejected. The triangle $(14,30,40)$ shows that the problem is not merely about right triangles. Its sines $7/25$, $3/5$, $4/5$ come from the primitive triples $7,24,25$ and $3,4,5$, and the angles do close up: with $A,B,C$ opposite $14,30,40$ the cosines are $24/25$, $4/5$, $-3/5$, so the addition formulas give $\sin(A+B)=\sin C$ and $\cos(A+B)=-\cos C$; as $A$ and $B$ are acute, $A+B$ lies in $(0,\pi)$, which forces $A+B=\pi-C$.

9. Which radii contribute at all

If $r$ has no prime factor congruent to $1$ modulo $4$ then $\mathcal S_r=\{2r\}$, and the only triple available is the equilateral $(2r,2r,2r)$, whose circumradius is $2r/\sqrt3$, never $r$: such radii contribute nothing.

Conversely, if some $d\mid r$ is a primitive hypotenuse, say $d=m^2+n^2$, then with $k=2r/d$ the three candidates $k(m^2-n^2)$, $2kmn$ and $2r=kd$ satisfy $(m^2-n^2)^2+(2mn)^2=d^2$ after dividing by $k^2$, so they form a right triangle whose hypotenuse is a diameter and whose circumradius is exactly $r$. Hence $T(r)\ge1$ precisely when $r$ has a prime factor congruent to $1$ modulo $4$. Of the radii up to $10^7$, exactly $7449349$ are of this kind. The code has no early exit for the other $2550651$: it still factors them and walks their divisors, but the candidate list is only $\{2r\}$, so a single triple reaches the exact test and fails.

10. Why the count is complete and has no repetitions

A triple that passes the triangle inequality and the identity of Section 1 with $R=r$ is a genuine triangle of circumradius exactly $r$; conversely Sections 2 to 6 place all three sides of any such triangle in $\mathcal S_r$; and the enumeration runs over $a\le b\le c$ in sorted order, so each multiset of sides — each congruence class — is offered exactly once. Hence $\sum_r r\,T(r)$ is exactly $S(n)$.

How the Code Works

The products in the test of Section 1 are large: with $a,b,c\le2\cdot10^7$ the value $a^2b^2c^2$ can reach $6.4\cdot10^{43}$, which needs $146$ bits and so overflows $128$. The C++ file therefore defines a small fixed-width type U192 holding three 64-bit limbs lo, mid, hi, with a single operation multiply_small that multiplies the whole value by one 64-bit factor, propagating carries through unsigned __int128 intermediates, ^[9] and an equality comparison. circumradius_equals(a, b, c, radius) builds both sides of the identity from U192::one() by six such multiplications each — $r,r$ and the four bracketed factors on one side, $a,a,b,b,c,c$ on the other — and compares them. Nothing is divided and nothing is rounded.

The constructor of Euler373Solver fills two tables. build_smallest_prime_factor is a sieve that stores in smallest_prime_factor_[v] the least prime dividing $v$: it initialises each entry to $v$ and, for every $p$ with $p^2\le N$ still equal to its own entry, overwrites the untouched multiples of $p$ from $p^2$ upwards. Because $p$ increases, the smallest prime wins, and every composite $v\le N$ is reached since its least prime factor $p$ satisfies $p^2\le v$. build_denominator_table realises Section 5: it loops over $m$ up to $\lfloor\sqrt N\rfloor$ and $n\lt m$, skips pairs of equal parity or with $\gcd(m,n)\ne1$, leaves the inner loop as soon as $m^2+n^2\gt N$ — legitimate because $m^2+n^2$ grows with $n$ — and appends $m^2-n^2$ and $2mn$ to denominator_to_numerators_[m*m+n*n]. Afterwards each list is sorted and deduplicated.

divisors_of peels the prime powers of $r$ off with the sieve table and expands the divisor list prime by prime, multiplying the divisors found so far by $p,p^2,\dots,p^{e}$. count_triangles_for_radius then builds $\mathcal S_r$ exactly as in Section 7: it starts the vector sides with the diameter 2 * radius, and for each divisor denom with a non-empty numerator list it appends scale * numer with scale = 2 * (radius / denom). The list is sorted and passed through std::unique. Three nested loops over indices $i\le j\le k$ then walk the sorted candidates; the inner loop uses if (ab <= c) break, which is correct because $c$ only grows, so once the triangle inequality fails for one $c$ it fails for all later ones. Surviving triples go to circumradius_equals. sum_up_to adds $r$ times the resulting count into an unsigned __int128 accumulator, and since the contributions are independent the order of accumulation cannot change the exact total.

Two self-checks are packaged with the program. brute_force_triangle_count_for_radius ignores all of the theory and tries every $a\le b\le c\le 2r$ with $c\le a+b-1$, applying the same exact test; run_validation_checkpoints compares it with count_triangles_for_radius for every $r\le60$ and then recomputes $S(100)$ and $S(1200)$, aborting with a message on standard error if either differs. Only then does main build the solver for $10^7$ and print the result, converted to decimal by to_string_u128 because the standard streams cannot format a 128-bit value.

The Python and Java files are wrappers, not independent implementations. Each expects solutionsCpp/Euler373.cpp under its project root — the parent of the script directory in Python, the working directory in Java — compiles it with clang++ or g++ at -std=c++17 -O2 into a cached bridge binary, recompiling only when the source is newer, and runs it. parse_output in Python and parseOutput in Java then take the last value after answer:, failing that the last after =, failing that the last non-empty line; neither checks that it is numeric, and here the last branch applies. The mathematics, the exact test and the checkpoints are thus shared by all three entry points. Both also look for --skip-checkpoints in the C++ source and would forward it; that string does not occur here and main takes no arguments, so the self-checks always run.

Complexity and Verification

Let $N=10^7$. The sieve costs $O(N\log\log N)$ time and one int per value. The parameter search of build_denominator_table visits $3926117$ pairs $(m,n)$, which is $O(N)$ pairs since both parameters stay below $\sqrt N$, each costing one gcd, so this stage is $O(N\log N)$. It keeps $1591579$ primitive triples, matching Lehmer's asymptotic count $N/(2\pi)\approx1591549$ for primitive triples with hypotenuse below $N$ to within $30$; ^[8] they are spread over $814182$ distinct hypotenuses and stored as $3183158$ numerators.

Writing $\tau(r)$ for the number of divisors of $r$, factoring a radius with the sieve table and expanding its divisors costs $O(\tau(r)+\log r)$, and $\sum_{r\le N}\tau(r)=162725364$, close to the classical estimate $N\ln N+(2\gamma-1)N\approx162725270$ with $\gamma\approx0.577216$ the Euler–Mascheroni constant. ^[10] Building $\mathcal S_r$ scans all $\tau(r)$ divisors, nearly all of which carry no numerator list, and then sorts the candidates, so one radius costs $O(\tau(r)+|\mathcal S_r|\log|\mathcal S_r|)$. The candidate total is $\sum_{r\le N}|\mathcal S_r|=56942950$, a mean of only $5.69$ per radius. The triple scan dominates everything else: it is $O(|\mathcal S_r|^3)$ in the worst case, and over all radii it performs $5453418759$ exact tests, each of which makes $12$ calls to multiply_small with three limb products apiece: about $2\cdot10^{11}$ limb multiplications at source level, whose instruction count depends on the compiler and the target. This is why $|\mathcal S_r|$ matters: the single radius $5928325$, with $405$ candidates, accounts for $7188057$ of those tests. Memory is dominated by the two tables, measured at about $300$ megabytes, of which the vector of $N+1$ numerator lists takes some $240$ in list headers alone.

Beyond the packaged checks — brute force against the fast enumeration for every $r\le60$, and $S(100)=4950$, $S(1200)=1653605$ — the following were verified for this article. An independent brute force agrees with the candidate enumeration for every $r\le120$, past the range the program tests, and the method further gives $S(10^4)=200174860$ and $S(10^5)=33211076049$. Over all $r\le10^7$ the deduplication of sides never removes an element, confirming the distinctness argument of Section 7, and $|\mathcal S_r|=\prod(1+2e)$ matches the constructed lists for every $r\le20000$. Up to $10^5$ the non-empty numerator lists sit exactly at the integers above $1$ whose prime factors are all congruent to $1$ modulo $4$, ^[7] every stored leg is coprime to its hypotenuse, every primitive hypotenuse found is odd, and every radius with more than one candidate does yield a triangle, as Section 9 predicts.

Footnotes and References

The references below supply the background facts used above. The reduction of the sides to legs of primitive triples whose hypotenuse divides the radius, the count $|\mathcal S_r|$, the characterisation of the contributing radii and the operation counts are derived or measured in this article.

Project Euler 373 — Circumscribed Circles. The official statement restricts attention to integer-sided triangles whose circumradius is integral as well, defines S(n) as the sum of the radii of all such triangles whose radius does not exceed n, gives S(100) = 4950 and S(1200) = 1653605, and asks for S(10⁷).
Wikipedia — Circumcircle. The page gives the diameter of a triangle's circumscribed circle as abc/(2·Area), and also as 2abc divided by the square root of (a+b+c)(−a+b+c)(a−b+c)(a+b−c). Halving these gives R = abc/(4Δ) as used in Section 1, and the same identity with Heron's product substituted.
Wikipedia — Heron's formula. The area is one quarter of the square root of (a+b+c)(−a+b+c)(a−b+c)(a+b−c), which squared is the form used in Section 1. It also states that lengths obeying the strict triangle inequality give a positive real area, that equality gives zero area and that violating it makes the formula imaginary — the sign facts behind the bracketed factors.
Wikipedia — Law of sines. In its extended form the common value of a/sin α, b/sin β and c/sin γ is 2R, the diameter of the circumcircle — the relation a = 2r sin A of Section 2. The page also gives the area as abc/(4R), the identity used in Section 1.
Wikipedia — Law of cosines. The relation c² = a² + b² − 2ab cos γ, solved for the cosine, makes the cosine of every angle a quotient of integers when the sides are integers, which is what Section 3 needs.
Wikipedia — Pythagorean triple. The page parameterises the primitive triples by coprime m > n > 0 of opposite parity as m² − n², 2mn and m² + n², with the two legs possibly exchanged — which is why the table of Section 5 stores both numerators. It also states that a primitive hypotenuse is always odd, that an integer greater than 1 is one if and only if each of its prime factors is congruent to 1 modulo 4, and that a hypotenuse with k distinct prime factors carries 2^(k−1) primitive triples.
OEIS A008846 — Hypotenuses of primitive Pythagorean triangles. The sequence begins 5, 13, 17, 25, 29, 37, 41, 53, 61, 65, lists these numbers as the values x² + y² with x even, y odd and gcd(x, y) = 1, and notes that all divisors of its terms are of the form 4k + 1.
MathWorld — Pythagorean Triple. The page reports Lehmer's 1900 result that the number of primitive triples with hypotenuse below N divided by N tends to 1/(2π) = 0.1591549…, and tabulates that count as 1, 16, 158 and 1593 for N = 10, 100, 1000 and 10000 — values matched by the table the program builds.
GCC manual — 128-bit integers. The page documents __int128 and unsigned __int128 as an extension available on targets with an integer mode wide enough for 128 bits. The three-limb type uses them to catch the carries of each 64-bit multiplication, and the final sum is accumulated in one.
Wikipedia — Divisor summatory function. The sum of the number of divisors over all n up to x equals x log x + (2γ − 1)x plus an error that Dirichlet bounded by O(√x).

Problem 373 source code

C++

#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <string>
#include <thread>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

constexpr int kTargetRadius = 10'000'000;

struct U192 {
    u64 lo = 0;
    u64 mid = 0;
    u64 hi = 0;

    static U192 one() {
        U192 value;
        value.lo = 1;
        return value;
    }

    void multiply_small(u64 factor) {
        const u128 t0 = static_cast<u128>(lo) * factor;
        lo = static_cast<u64>(t0);
        u64 carry = static_cast<u64>(t0 >> 64);

        const u128 t1 = static_cast<u128>(mid) * factor + carry;
        mid = static_cast<u64>(t1);
        carry = static_cast<u64>(t1 >> 64);

        const u128 t2 = static_cast<u128>(hi) * factor + carry;
        hi = static_cast<u64>(t2);
    }

    bool operator==(const U192& other) const {
        return lo == other.lo && mid == other.mid && hi == other.hi;
    }
};

std::string to_string_u128(u128 value) {
    if (value == 0) {
        return "0";
    }

    std::string digits;
    while (value > 0) {
        digits.push_back(static_cast<char>('0' + value % 10));
        value /= 10;
    }

    std::reverse(digits.begin(), digits.end());
    return digits;
}

bool circumradius_equals(int a, int b, int c, int radius) {
    // Exact check of:
    // a^2 b^2 c^2 == R^2 (a+b+c)(a+b-c)(a-b+c)(-a+b+c)
    U192 lhs = U192::one();
    lhs.multiply_small(static_cast<u64>(radius));
    lhs.multiply_small(static_cast<u64>(radius));
    lhs.multiply_small(static_cast<u64>(a + b + c));
    lhs.multiply_small(static_cast<u64>(a + b - c));
    lhs.multiply_small(static_cast<u64>(a - b + c));
    lhs.multiply_small(static_cast<u64>(-a + b + c));

    U192 rhs = U192::one();
    rhs.multiply_small(static_cast<u64>(a));
    rhs.multiply_small(static_cast<u64>(a));
    rhs.multiply_small(static_cast<u64>(b));
    rhs.multiply_small(static_cast<u64>(b));
    rhs.multiply_small(static_cast<u64>(c));
    rhs.multiply_small(static_cast<u64>(c));

    return lhs == rhs;
}

class Euler373Solver {
public:
    explicit Euler373Solver(int radius_limit)
        : radius_limit_(radius_limit),
          smallest_prime_factor_(static_cast<std::size_t>(radius_limit) + 1),
          denominator_to_numerators_(static_cast<std::size_t>(radius_limit) + 1) {
        build_smallest_prime_factor();
        build_denominator_table();
    }

    u128 solve(bool allow_multithreading, unsigned requested_threads = 0) const {
        return sum_up_to(radius_limit_, allow_multithreading, requested_threads);
    }

    u128 sum_up_to(int max_radius,
                   bool allow_multithreading,
                   unsigned requested_threads = 0) const {
        if (max_radius > radius_limit_) {
            return 0;
        }

        unsigned threads = requested_threads;
        if (threads == 0) {
            threads = std::thread::hardware_concurrency();
            if (threads == 0) {
                threads = 1;
            }
        }

        if (!allow_multithreading || threads <= 1 || max_radius < 50'000) {
            threads = 1;
        } else {
            threads = std::min<unsigned>(threads, static_cast<unsigned>(max_radius));
        }

        if (threads == 1) {
            std::vector<int> divisors;
            std::vector<int> sides;
            u128 total = 0;
            for (int radius = 1; radius <= max_radius; ++radius) {
                total += static_cast<u128>(radius) *
                         count_triangles_for_radius(radius, divisors, sides);
            }
            return total;
        }

        std::atomic<int> next_radius{1};
        std::vector<u128> partial(threads, 0);
        std::vector<std::thread> workers;
        workers.reserve(threads);

        for (unsigned tid = 0; tid < threads; ++tid) {
            workers.emplace_back([&, tid]() {
                std::vector<int> divisors;
                std::vector<int> sides;
                u128 local = 0;

                while (true) {
                    const int radius = next_radius.fetch_add(1, std::memory_order_relaxed);
                    if (radius > max_radius) {
                        break;
                    }

                    local += static_cast<u128>(radius) *
                             count_triangles_for_radius(radius, divisors, sides);
                }

                partial[tid] = local;
            });
        }

        for (std::thread& worker : workers) {
            worker.join();
        }

        u128 total = 0;
        for (const u128 value : partial) {
            total += value;
        }
        return total;
    }

    u64 count_triangles_for_radius(int radius) const {
        std::vector<int> divisors;
        std::vector<int> sides;
        return count_triangles_for_radius(radius, divisors, sides);
    }

private:
    int radius_limit_ = 0;
    std::vector<int> smallest_prime_factor_;
    std::vector<std::vector<int>> denominator_to_numerators_;

    void build_smallest_prime_factor() {
        for (int value = 0; value <= radius_limit_; ++value) {
            smallest_prime_factor_[static_cast<std::size_t>(value)] = value;
        }
        if (radius_limit_ >= 1) {
            smallest_prime_factor_[1] = 1;
        }

        for (int p = 2; static_cast<std::int64_t>(p) * p <= radius_limit_; ++p) {
            if (smallest_prime_factor_[static_cast<std::size_t>(p)] != p) {
                continue;
            }

            for (int multiple = p * p; multiple <= radius_limit_; multiple += p) {
                if (smallest_prime_factor_[static_cast<std::size_t>(multiple)] == multiple) {
                    smallest_prime_factor_[static_cast<std::size_t>(multiple)] = p;
                }
            }
        }
    }

    void build_denominator_table() {
        const int m_limit = static_cast<int>(std::sqrt(static_cast<long double>(radius_limit_)));

        for (int m = 2; m <= m_limit; ++m) {
            for (int n = 1; n < m; ++n) {
                if (((m + n) & 1) == 0) {
                    continue;
                }
                if (std::gcd(m, n) != 1) {
                    continue;
                }

                const int denom = m * m + n * n;
                if (denom > radius_limit_) {
                    break;
                }

                auto& numerators = denominator_to_numerators_[static_cast<std::size_t>(denom)];
                numerators.push_back(m * m - n * n);
                numerators.push_back(2 * m * n);
            }
        }

        for (int denom = 1; denom <= radius_limit_; ++denom) {
            auto& numerators = denominator_to_numerators_[static_cast<std::size_t>(denom)];
            if (numerators.empty()) {
                continue;
            }
            std::sort(numerators.begin(), numerators.end());
            numerators.erase(std::unique(numerators.begin(), numerators.end()),
                             numerators.end());
        }
    }

    void divisors_of(int value, std::vector<int>& out) const {
        out.clear();
        out.push_back(1);

        while (value > 1) {
            const int prime = smallest_prime_factor_[static_cast<std::size_t>(value)];
            int exponent = 0;
            while (value % prime == 0) {
                value /= prime;
                ++exponent;
            }

            const std::size_t base_size = out.size();
            int multiplier = 1;
            for (int i = 1; i <= exponent; ++i) {
                multiplier *= prime;
                for (std::size_t idx = 0; idx < base_size; ++idx) {
                    out.push_back(out[idx] * multiplier);
                }
            }
        }
    }

    u64 count_triangles_for_radius(int radius,
                                   std::vector<int>& divisors,
                                   std::vector<int>& sides) const {
        divisors_of(radius, divisors);

        sides.clear();
        sides.push_back(2 * radius); // Diameter (sin(theta)=1).

        for (const int denom : divisors) {
            const auto& numerators =
                denominator_to_numerators_[static_cast<std::size_t>(denom)];
            if (numerators.empty()) {
                continue;
            }

            const int scale = 2 * (radius / denom);
            for (const int numer : numerators) {
                sides.push_back(scale * numer);
            }
        }

        std::sort(sides.begin(), sides.end());
        sides.erase(std::unique(sides.begin(), sides.end()), sides.end());

        const int count = static_cast<int>(sides.size());
        u64 triangles = 0;

        for (int i = 0; i < count; ++i) {
            const int a = sides[static_cast<std::size_t>(i)];
            for (int j = i; j < count; ++j) {
                const int b = sides[static_cast<std::size_t>(j)];
                const int ab = a + b;
                for (int k = j; k < count; ++k) {
                    const int c = sides[static_cast<std::size_t>(k)];
                    if (ab <= c) {
                        break;
                    }

                    if (circumradius_equals(a, b, c, radius)) {
                        ++triangles;
                    }
                }
            }
        }

        return triangles;
    }
};

u64 brute_force_triangle_count_for_radius(int radius) {
    const int max_side = 2 * radius;
    u64 count = 0;

    for (int a = 1; a <= max_side; ++a) {
        for (int b = a; b <= max_side; ++b) {
            const int c_max = std::min(max_side, a + b - 1);
            for (int c = b; c <= c_max; ++c) {
                if (circumradius_equals(a, b, c, radius)) {
                    ++count;
                }
            }
        }
    }

    return count;
}

bool run_validation_checkpoints() {
    {
        // Structural validation: compare the optimized candidate enumeration
        // against brute force on small radii.
        Euler373Solver solver(60);
        for (int radius = 1; radius <= 60; ++radius) {
            const u64 fast = solver.count_triangles_for_radius(radius);
            const u64 brute = brute_force_triangle_count_for_radius(radius);
            if (fast != brute) {
                std::cerr << "Validation failed for radius " << radius << ": fast=" << fast
                          << ", brute=" << brute << '\n';
                return false;
            }
        }
    }

    struct Checkpoint {
        int radius_limit;
        u64 expected;
    };

    const std::vector<Checkpoint> checkpoints = {
        {100, 4'950},
        {1200, 1'653'605},
    };

    for (const Checkpoint cp : checkpoints) {
        const Euler373Solver solver(cp.radius_limit);
        const u64 got = static_cast<u64>(solver.solve(false, 1));
        if (got != cp.expected) {
            std::cerr << "Checkpoint failed for S(" << cp.radius_limit << "): got " << got
                      << ", expected " << cp.expected << '\n';
            return false;
        }
    }

    return true;
}

} // namespace

int main() {
    if (!run_validation_checkpoints()) {
        return 1;
    }

    unsigned threads = std::thread::hardware_concurrency();
    if (threads == 0) {
        threads = 4;
    }

    const Euler373Solver solver(kTargetRadius);
    const u128 answer = solver.solve(true, threads);
    std::cout << to_string_u128(answer) << '\n';
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

import java.nio.file.*;
import java.util.*;
import java.util.regex.*;

public class Euler373 {
    private static final Pattern ANSWER_RE = Pattern.compile("answer\\s*:\\s*(.+)$", Pattern.CASE_INSENSITIVE);
    private static final Pattern EQUAL_RE = Pattern.compile("=\\s*(.+)$");

    private static String parseOutput(String stdout) {
        String[] lines = stdout.split("\\R");
        List<String> nonEmpty = new ArrayList<>();
        for (String line : lines) {
            String t = line.trim();
            if (!t.isEmpty()) {
                nonEmpty.add(t);
            }
        }
        if (nonEmpty.isEmpty()) {
            return "";
        }

        List<String> answers = new ArrayList<>();
        List<String> equals = new ArrayList<>();
        for (String line : nonEmpty) {
            Matcher m1 = ANSWER_RE.matcher(line);
            if (m1.find()) {
                answers.add(m1.group(1).trim());
            }
            Matcher m2 = EQUAL_RE.matcher(line);
            if (m2.find()) {
                equals.add(m2.group(1).trim());
            }
        }

        if (!answers.isEmpty()) {
            return answers.get(answers.size() - 1);
        }
        if (!equals.isEmpty()) {
            return equals.get(equals.size() - 1);
        }
        return nonEmpty.get(nonEmpty.size() - 1);
    }

    private static String pickCompiler() throws Exception {
        for (String compiler : List.of("clang++", "g++")) {
            Process probe = new ProcessBuilder("bash", "-lc", "command -v " + compiler)
                    .redirectErrorStream(true)
                    .start();
            String out = new String(probe.getInputStream().readAllBytes());
            int rc = probe.waitFor();
            if (rc == 0 && !out.trim().isEmpty()) {
                return compiler;
            }
        }
        throw new RuntimeException("No C++ compiler found (clang++/g++).");
    }

    private static Path cppSource(Path root) {
        return root.resolve("solutionsCpp").resolve("Euler373.cpp");
    }

    private static boolean shouldSkipCheckpoints(Path root) {
        Path src = cppSource(root);
        try {
            String text = Files.readString(src);
            return text.contains("--skip-checkpoints");
        } catch (Exception ex) {
            return false;
        }
    }

    private static Path ensureBridgeBinary() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = root.resolve("solutionsCpp").resolve(".euler373_java_bridge");

        boolean rebuild = Files.notExists(bin)
                || Files.getLastModifiedTime(src).compareTo(Files.getLastModifiedTime(bin)) > 0;

        if (rebuild) {
            String compiler = pickCompiler();
            Process compile = new ProcessBuilder(
                    compiler,
                    "-std=c++17",
                    "-O2",
                    src.toString(),
                    "-o",
                    bin.toString())
                    .inheritIO()
                    .start();
            if (compile.waitFor() != 0) {
                throw new RuntimeException("Failed to compile Euler373 C++ bridge.");
            }
        }

        return bin;
    }

    private static String runBridge(Path bin, Path root, Path srcDir) throws Exception {
        List<String> cmd = new ArrayList<>();
        cmd.add(bin.toString());
        if (shouldSkipCheckpoints(root)) {
            cmd.add("--skip-checkpoints");
        }

        Process first = new ProcessBuilder(cmd)
                .directory(root.toFile())
                .redirectErrorStream(true)
                .start();
        String out = new String(first.getInputStream().readAllBytes());
        int rc = first.waitFor();
        if (rc == 0) {
            return out;
        }

        Process second = new ProcessBuilder(cmd)
                .directory(srcDir.toFile())
                .redirectErrorStream(true)
                .start();
        String out2 = new String(second.getInputStream().readAllBytes());
        int rc2 = second.waitFor();
        if (rc2 == 0) {
            return out2;
        }

        throw new RuntimeException("Euler373 C++ bridge failed.\n" + out + "\n" + out2);
    }

    private static String solveViaCppBridge() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = ensureBridgeBinary();
        String out = runBridge(bin, root, src.getParent());
        String parsed = parseOutput(out);
        if (parsed.isEmpty()) {
            throw new RuntimeException("Euler373 C++ bridge produced empty output.");
        }
        return parsed;
    }

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