Problem 143: Investigating the Torricelli Point of a Triangle

Project Euler Problem 143 Solution

EulerSolve provides an optimized solution for Project Euler Problem 143, Investigating the Torricelli Point of a Triangle, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $T$ be the Torricelli point of a triangle $ABC$, and write $p=TA$, $q=TB$, $r=TC$. The problem asks for the sum of all distinct values of $p+q+r$ with $p+q+r \le 120{,}000$, under the condition that $p$, $q$, and $r$ are integers and that the three side lengths of the triangle are integers as well. The decisive fact is that the three angles at $T$ are all $120^\circ$. That turns the geometry into an arithmetic question about when expressions of the form $x^2+xy+y^2$ are perfect squares. The implementations exploit that arithmetic structure directly and never need to reconstruct coordinates. Mathematical Approach Write the triangle sides as $a=BC$, $b=CA$, and $c=AB$. Since the angles $\angle ATB$, $\angle BTC$, and $\angle CTA$ are each $120^\circ$, every side length is determined by a cosine-law identity. From the Torricelli point to a Diophantine system Applying the law of cosines with $\cos 120^\circ=-\tfrac12$ gives $$a^2=q^2+qr+r^2,\qquad b^2=r^2+rp+p^2,\qquad c^2=p^2+pq+q^2.$$ So a triple $(p,q,r)$ is admissible exactly when all three quadratic forms $$p^2+pq+q^2,\qquad q^2+qr+r^2,\qquad r^2+rp+p^2$$ are perfect squares. Once this happens, the corresponding triangle sides are integers, and the point with three $120^\circ$ rays is precisely the Torricelli configuration described in the problem....

Detailed mathematical approach

Problem Summary

Let $ABC$ be a triangle whose angles are all less than $120^\circ$, with sides $a=BC$, $b=CA$ and $c=AB$, and for a point $X$ inside it put $p=XA$, $q=XC$ and $r=XB$. Fermat challenged Torricelli to find the point $X$ that minimises $p+q+r$. Torricelli proved that the circumcircles of the equilateral triangles $AOB$, $BNC$ and $AMC$ constructed on the sides meet in one point $T$ inside the triangle, the Torricelli point, and that $T$ minimises the sum. The triangle is a Torricelli triangle if, with the sum minimised, $a,b,c,p,q,r$ are all positive integers, as for $a=399$, $b=455$, $c=511$ with $p+q+r=784$. We must find the sum of all distinct values of $p+q+r\le120000$ over all Torricelli triangles. ^[1]

The whole solution in one sentence: from $T$ every side is seen under an angle of $120^\circ$, so by the law of cosines the Torricelli triangles correspond to the triples of positive integers $p,q,r$ for which $p^2+pq+q^2$, $q^2+qr+r^2$ and $r^2+rp+p^2$ are all squares; a two-parameter formula and its multiples list every pair of numbers with this property exactly once, and we add up the distinct sums of the triples whose three pairs all occur in that list.

Mathematical Approach

We turn the geometry into three conditions on whole numbers, find a formula for all pairs that satisfy one of them, and search for the triples in which all three pairs do.

1. Three angles of $120^\circ$

The equilateral triangles of the construction stand outside $ABC$, so the apex $O$ on $AB$ and the point $T$ inside $ABC$ lie on opposite sides of the line $AB$. As $T$ lies on the circle through $A$, $O$ and $B$, the quadrilateral $AOBT$ is inscribed in a circle, and its opposite angles add up to $180^\circ$. ^[3] So $\angle ATB=180^\circ-\angle AOB$${}=120^\circ$, and likewise $\angle BTC=\angle CTA=120^\circ$. ^[2]

Claim. If every side of a triangle is seen under $120^\circ$ from a point $T$ inside it, then $XA+XB+XC\gt TA+TB+TC$ for every point $X\ne T$.

Proof. Let $u_A$, $u_B$, $u_C$ be the unit vectors from $T$ towards $A$, $B$, $C$. A rotation by $120^\circ$ about $T$ permutes them, so it does not change their sum, and only the zero vector is unchanged by such a rotation: $u_A+u_B+u_C=0$. The dot product of a unit vector with a vector is at most the length of that vector, and $A-T=TA\cdot u_A$, so for every point $X$

$$\begin{aligned}XA&=|A-X|\\&\ge u_A\cdot(A-X)\\&=TA+u_A\cdot(T-X).\end{aligned}$$

When the three inequalities of this kind are added, the terms with $T-X$ cancel because $u_A+u_B+u_C=0$, and $XA+XB+XC\ge TA+TB+TC$ remains. Equality requires $A-X$ to be a nonnegative multiple of $u_A$, so $X$ lies on the line $TA$, and likewise on the lines $TB$ and $TC$, which meet only at $T$. ^[2] So the minimising point is unique, and in a Torricelli triangle it is $T$.

2. From the triangle to three squares

The law of cosines in the triangle $ATB$ gives $c^2=p^2+r^2-2pr\cos120^\circ$, ^[4] and $\cos120^\circ=-\tfrac12$. With the triangles $BTC$ and $CTA$ we get

$$\begin{aligned}a^2&=q^2+qr+r^2,\\b^2&=p^2+pq+q^2,\\c^2&=p^2+pr+r^2.\end{aligned}$$

In the example $p=325$, $q=195$ and $r=264$: $195^2+195\cdot264+264^2$${}=159201=399^2$, $325^2+325\cdot195+195^2$${}=207025=455^2$ and $325^2+325\cdot264+264^2$${}=261121=511^2$.

Conversely, let $p,q,r$ be positive integers for which the three right-hand sides are squares $a^2,b^2,c^2$. Draw three rays from a point $T$ at angles of $120^\circ$ to each other, and mark $A$, $C$, $B$ on them at the distances $p$, $q$, $r$. With $T$ as the origin, $u_A+u_B+u_C=0$ reads $A/p+B/r+C/q=0$, so $T$ is an average of the vertices with positive weights and lies inside $ABC$. The triangle $ATB$ has the angle $120^\circ$ at $T$, so its other two angles are below $60^\circ$; hence the angle $\angle TAB+\angle TAC$ of $ABC$ at $A$ is below $120^\circ$, and so are the other two. The sides are $a$, $b$, $c$, and by the claim $T$ minimises the sum, so $ABC$ is a Torricelli triangle. We therefore need the distinct sums $p+q+r\le120000$ of such triples.

3. Compatible pairs

Call positive integers $x$ and $y$ compatible if $x^2+xy+y^2=z^2$ for an integer $z$, the third side of a triangle in which $x$ and $y$ enclose an angle of $120^\circ$. A Torricelli triangle is the same as three pairwise compatible numbers. The compatible pairs with both numbers at most $40$ are

$(3,5)$, $(5,16)$, $(6,10)$, $(7,8)$, $(7,33)$, $(9,15)$, $(10,32)$, $(11,24)$, $(12,20)$, $(13,35)$, $(14,16)$, $(15,25)$, $(16,39)$, $(18,30)$, $(21,24)$, $(21,35)$, $(24,40)$, $(28,32)$, $(35,40)$.

Seven of them have no common factor, $(3,5)$, $(7,8)$, $(5,16)$, $(11,24)$, $(7,33)$, $(13,35)$ and $(16,39)$, with $z=7,13,19,31,37,43,49$; the others are multiples, such as $(6,10)=2\cdot(3,5)$. Two facts follow from the definition.

No number is compatible with itself: $3x^2$ contains the prime $3$ an odd number of times, every square an even number of times. So $p$, $q$, $r$ are different, and since the three conditions treat them alike, from now on we sort them as $p\lt q\lt r$.
Common factors can be removed: $(kx,ky)$ is compatible with the third side $kz$; conversely, for $g=\gcd(x,y)$, $g^2\mid z^2$ gives $g\mid z$, so $(x/g,y/g)$ is a compatible pair without a common factor, a primitive pair.

4. A formula that produces compatible pairs

For integers $m\gt n\ge1$ put

$$\begin{aligned}x_0&=m^2-n^2,\\y_0&=n(2m+n),\\z_0&=m^2+mn+n^2.\end{aligned}$$

Multiplying out gives $x_0^2+x_0y_0+y_0^2=z_0^2$. The reason: the Eisenstein integers $u+v\omega$, with integers $u$, $v$ and $\omega=(-1+i\sqrt3)/2$, satisfy $|u+v\omega|^2=u^2-uv+v^2$, ^[6] so $x^2+xy+y^2=|x-y\omega|^2$; and since $\omega^2=-1-\omega$,

$$(m-n\omega)^2=x_0-y_0\,\omega,$$

so $|x_0-y_0\omega|=|m-n\omega|^2=z_0$, and squaring gives the identity. The pairs $(m,n)=(2,1)$, $(3,1)$, $(3,2)$, $(5,1)$ give $(x_0,y_0,z_0)=(3,5,7)$, $(8,7,13)$, $(5,16,19)$, $(24,11,31)$. The pair $(4,1)$ gives $(15,9,21)=3\cdot(5,3,7)$, the first triple again, exchanged and tripled; for coprime $m$ and $n$, the next section shows that this happens exactly when $3$ divides $m-n$. The same formula describes the integer triangles with an angle of $120^\circ$. ^[5]

5. Every compatible pair arises exactly once

Call $(m,n)$ a seed if $m\gt n\ge1$, $\gcd(m,n)=1$ and $3\nmid m-n$.

Claim. Every compatible pair, in one of its two orders, equals $(kx_0,ky_0)$ for exactly one seed $(m,n)$ and one $k\ge1$.

Proof. Let $x^2+xy+y^2=z^2$ with $x,y\ge1$, and put $t=y/(x+z)$; since $z\gt y$, we have $0\lt t\lt1$. From $(z-x)(z+x)=y(x+y)$ we get $z-x=t(x+y)$, and $y=t(z+x)$ by definition. Solving these two linear equations gives

$$\begin{aligned}\frac xz&=\frac{1-t^2}{1+t+t^2},\\ \frac yz&=\frac{t(2+t)}{1+t+t^2},\end{aligned}$$

so with $t=n/m$ in lowest terms, $(x,y,z)$ is a positive rational multiple of $(x_0,y_0,z_0)$. Here $t$ is the slope of the line from $(-1,0)$ to the point $(x/z,y/z)$ of the curve $X^2+XY+Y^2=1$: lines through one known rational point find all rational points of such a curve. ^[7]

A prime $d$ dividing $x_0$, $y_0$ and $z_0$ cannot divide $n$, as it would then divide $m^2=x_0+n^2$; so $d\mid2m+n$. It divides one of the factors $m-n$ and $m+n$ of $x_0$, but not $m+n$, as it would then divide $(2m+n)-(m+n)=m$. So $d$ divides $(2m+n)+(m-n)=3m$ but not $m$; hence $d=3$ and $3\mid m-n$. For a seed, $x_0$, $y_0$, $z_0$ therefore have no common factor. If the rational multiple is $e/f$ in lowest terms, $f$ divides $ex_0$, $ey_0$ and $ez_0$, hence $x_0$, $y_0$ and $z_0$, so $f=1$: the multiple is an integer $k$, and $(x,y)=(kx_0,ky_0)$.

Exchanging $x$ and $y$ turns $t$ into $x/(y+z)$, and by the formulas above

$$\begin{aligned}\frac{x}{y+z}&=\frac{m^2-n^2}{m^2+3mn+2n^2}\\&=\frac{m-n}{m+2n}.\end{aligned}$$

A common divisor of $m-n$ and $m+2n$ divides $3n$ and $3m$, so it is $1$ or $3$, and it is $3$ exactly when $3\mid m-n$. In that case the exchanged order has, in lowest terms, the parameters $((m+2n)/3,(m-n)/3)$, whose difference $n$ is not divisible by $3$, since otherwise $3$ would divide $m$ too: a seed. Otherwise it has the parameters $(m+2n,m-n)$, with the difference $3n$: not a seed. So exactly one order comes from a seed, and since scaling does not change $t$, that seed is unique, and so is $k$. For example, $(5,3)$ has $t=3/12=1/4$, leading to $(4,1)$, which is not a seed, while $(3,5)$ has $t=5/10=1/2$ and the seed $(2,1)$.

6. The example taken apart

Sorted, the distances of the example are $195\lt264\lt325$. The pair $(325,264)$ has $t=264/(325+511)=6/19$ and comes from the seed $(19,6)$ with $k=1$: $19^2-6^2=325$, $6\cdot(38+6)=264$ and $19^2+19\cdot6+6^2=511$. The pair $(195,264)$ is $3\cdot(65,88)$, and $(65,88)$, with $t=88/(65+133)=4/9$, comes from the seed $(9,4)$; so $k=3$ and $z=3\cdot133=399$. The pair $(195,325)$ is $65\cdot(3,5)$, from the seed $(2,1)$ with $k=65$ and $z=65\cdot7=455$.

7. Triangles in a graph

Fix a limit $L$, here $120000$, and join every two compatible numbers among $1,\dots,L$ by an edge. Three pairwise compatible numbers form a triangle of this graph, a clique of three vertices. ^[8] In a triangle $p\lt q\lt r$ with $p+q+r\le L$, every two of the numbers add up to less than $L$, and $r$ is a common neighbour of $p$ and $q$ that is larger than $q$. With every neighbour list sorted, the common neighbours of two vertices are found by walking along both lists together: compare the current entries, move on in the list with the smaller one, and record a match when they are equal. This takes at most as many steps as the two lists have entries.

Worked example, $L=1000$: the neighbours of $195$ are $29$, $80$, $117$, $264$, $325$, $525$, $624$, $952$, and those of $264$ are $56$, $121$, $195$, $231$, $325$, $440$, $576$, $945$. For the edge $195\lt264$ the walk starts after $264$, with $325,525,624,952$ and $325,440,576,945$. It finds $325$, the triangle of the example with the sum $784$, and then compares $525$ with $440$ and $576$, $624$ with $576$ and $945$, and $952$ with $945$, which ends the second list. The whole search for $L=1000$ finds the triangles $\{195,264,325\}$, $\{264,325,440\}$ and $\{390,528,650\}$, but only the first sum is at most $1000$; the others are $1029$ and $1568$. So $784$ is the only value up to $1000$, and the values up to $1500$ are $784$ and $1029$, with the total $1813$.

8. Distinct sums

Multiplying a Torricelli triple by $k$ multiplies its sum by $k$, so the values come in families. Those up to $3200$ are $784$, $1029$, $1568$, $2058$, $2352$, $2709$, $3087$ and $3136$: multiples of $784$ and $1029$, and $2709$ from the primitive triple $(384,805,1520)$. Different triangles can share a sum: $21\cdot(195,264,325)$ and $16\cdot(264,325,440)$ both have the sum $16464$, the least common multiple of $784=2^4\cdot7^2$ and $1029=3\cdot7^3$. As the problem asks for distinct values, $16464$ counts once: the total for $L=16464$ exceeds that for $L=16463$ by $16464$, not by twice as much. Up to $120000$ the $508$ triangles have only $491$ distinct sums, so the sums are collected in a set before they are added.

How the Code Works

The C++ program ^[10] starts in main. parse_arguments fills an Options structure, by default limit = 120000 and run_checkpoints = true; --limit=N is read by parse_int_after_prefix, which accepts only decimal digits and accumulates them in an int without an overflow check, and --skip-checkpoints turns the checks off. An unknown argument or a limit below $3$ gives exit status $1$, and a failed check in run_checkpoints status $2$. Otherwise main prints solve(limit), the sum of the set returned by generate_distinct_sums, added up in an unsigned 64-bit integer.

generate_distinct_sums stores the graph as a vector adj of $L+1$ lists of int. It runs over $2\le m\le\lfloor\sqrt{2L}\rfloor+4$ and $1\le n\lt m$, keeps only seeds (std::gcd(m, n) equal to $1$ and (m - n) % 3 not $0$), skips a seed whose x0 or y0 exceeds the limit, and for $k=1,2,\dots$ appends $kx_0$ and $ky_0$ to each other's lists while both are within the limit. The bound on $m$ is safe, since a seed inside the box has $m^2\lt x_0+y_0\le2L$. Each list is sorted and passed through std::unique, which by Section 5 finds nothing to remove. The search then takes every $p$ with at least two neighbours and every neighbour $q\gt p$ with $p+q\lt L$, finds the first entries above $q$ in both lists with std::upper_bound, ^[9] and walks the two tails as in Section 7; each common entry $r$ with $p+q+r\le L$ adds its sum to the std::set<int> sums. The walk does not stop early when the sums pass the limit. brute_small tests every $p\lt q$ with $p+q\le L$ and, for compatible ones, every $r\gt q$ with $p+q+r\le L$.

The Python version follows the C++ program function by function; it removes repeated entries with dict.fromkeys, finds the tails by a linear scan instead of a binary search, and uses integers that cannot overflow. The Java version has a fixed limit and no options or checks. It tests every pair $p\lt q$ with $p+q\le120000$, forming $n=p^2+pq+q^2$ with casts to long before the multiplications, since $pq$ can exceed the int range. isSquare compares the square of (long) Math.sqrt(n) with n; this is exact, because $n\lt(p+q)^2\le L^2\lt2^{53}$ is stored without rounding and Math.sqrt rounds correctly, ^[9] so the root of a perfect square comes out exactly. The pairs go into a HashMap of TreeSet objects. For every $p$ the program runs over the pairs of neighbours $q\lt r$, smaller ones than $p$ included, with the cut-offs $p+q\lt L$ and $p+q+r\le L$, and checks with contains whether $r$ is a neighbour of $q$; every triangle is found once from each corner, the TreeSet of sums keeps one copy, and the total is a long.

Complexity and Verification

For $L=120000$ the loop looks at $121278$ pairs $(m,n)$, about $L$ since $m_{\max}^2\approx2L$, each with a gcd; $21000$ of them are seeds inside the box, the largest with $m=371$. The larger of $x_0$, $y_0$ is at least their average, so a seed has at most $2L/(x_0+y_0)\lt2L/m^2$ multiples within the limit, the seeds with one $m$ give fewer than $2L/m$ pairs, and the graph has fewer than $2L\ln m_{\max}$ edges, which is $O(L\log L)$; in fact it has $220232$. Sorting costs $O(E\log D)$ for $E$ edges and the largest degree $D=126$. The search takes $140583$ edges $p\lt q$, each with two binary searches and a walk of at most $\deg p+\deg q$ steps, so the walks take at most $\sum_v(\deg v)^2=7\,304\,598$ steps; we know no closed formula for this sum. They actually take $1\,893\,497$ steps and find $886$ triangles, $508$ of them within the limit, each inserted into the set in $O(\log L)$ time. Python's linear scans stay within the same bound, with about $3.1$ million steps. The memory is $O(L+E)$. Counting operations on machine integers as single steps, this is a few million steps; on a desktop computer the C++ program needs about $0.03$ seconds with its checks, the Python program about $0.5$ seconds. The Java program makes $3\,599\,940\,000$ square tests, about $L^2/4$, which dwarf its $1\,654\,171$ lookups in the neighbour sets; it needs about $3$ seconds.

Before the main run, the C++ and Python programs check that $784$ is among the sums for the limit $1000$, that the total for $1000$ is exactly $784$, and that the total for $1500$ equals that of brute_small, which tests $561750$ pairs and finds the sums $784$ and $1029$. The Java program has no checks. Further checks, not part of the published solutions, were run for this article: a direct test of all $7\,199\,940\,000$ pairs $x\lt y\le120000$ found exactly the $220232$ pairs produced from the seeds, each once, as Section 5 predicts, and a triangle search on these pairs with hash lookups instead of sorted lists gave the same total as the published C++ program for every limit from $3$ to $3000$ and for $5000$, $16463$, $16464$, $20000$, $50000$ and $120000$. The C++ program, built with Clang and with GCC, the Python program and the Java program print the same value. The identities of Sections 4 and 5 were checked symbolically, and the claims about common factors and exchanged orders for all coprime $m\gt n$ with $m\lt400$.

Footnotes and References

The references below supply the background facts used above. The reduction to three square conditions, the proof that the seeds produce every compatible pair exactly once and the operation counts are derived in this article.

Project Euler 143 — Torricelli Triangles. The statement defines Torricelli triangles through the point T that minimises p + q + r, gives the example a = 399, b = 455, c = 511 with p + q + r = 784, and asks for the sum of all distinct values of p + q + r ≤ 120000.
Wikipedia — Fermat point. Fermat posed the problem to Torricelli, who solved it. If every angle of the triangle is less than 120°, the sides subtend 120° at the Fermat point and the circumcircles of the equilateral triangles on the sides pass through it; a vector proof of its minimality uses unit vectors with sum zero.
Wikipedia — Cyclic quadrilateral. A convex quadrilateral can be inscribed in a circle exactly when its opposite angles add up to 180°.
Wikipedia — Law of cosines. c² = a² + b² − 2ab cos γ, where γ is the angle opposite the side c.
Wikipedia — Integer triangle. Integer triangles with a 120° angle arise from a = m² + mn + n², b = 2mn + n², c = m² − n² with coprime 0 < n < m; gcd(a, b, c) is 3 if m ≡ n (mod 3) and 1 otherwise, and adding the condition m ≢ n (mod 3) generates every primitive triple once, a remark the page marks as needing a citation; Section 5 proves it.
Wikipedia — Eisenstein integer. The numbers a + bω with integers a, b and the cube root of unity ω = (−1 + i√3)/2 have the norm a² − ab + b², the square of their absolute value.
Wikipedia — Diophantine equation. If a non-trivial solution of a homogeneous quadratic equation is known and its point is not singular, the rational lines through it give all rational points; the article derives the Pythagorean triples this way from the point (−1, 0).
Wikipedia — Clique problem. A clique is a set of pairwise adjacent vertices; a graph with m edges has at most of the order of m^(3/2) triangles, and Chiba and Nishizeki list them all in O(m^(3/2)) time.
cppreference — std::upper_bound returns the first element ordered after the given value, with at most log₂ N + O(1) comparisons for N elements; Java SE 17 — Math.sqrt returns the correctly rounded positive square root of a double.
C++ — Euler143.cpp. The linked immutable revision contains generate_distinct_sums, solve, brute_small and run_checkpoints; it is the implementation record for this article.

Problem 143 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <set>
#include <string>
#include <vector>

namespace {

using u64 = std::uint64_t;

struct Options {
    int limit = 120000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return options.limit >= 3;
}

std::set<int> generate_distinct_sums(const int limit) {
    std::vector<std::vector<int>> adj(static_cast<std::size_t>(limit + 1));

    const int m_max = static_cast<int>(std::sqrt(2.0L * limit)) + 4;
    for (int m = 2; m <= m_max; ++m) {
        for (int n = 1; n < m; ++n) {
            if (std::gcd(m, n) != 1) {
                continue;
            }
            if ((m - n) % 3 == 0) {
                continue;
            }

            const int x0 = m * m - n * n;
            const int y0 = n * (2 * m + n);
            if (x0 <= 0 || y0 <= 0) {
                continue;
            }

            if (x0 > limit || y0 > limit) {
                continue;
            }

            for (int k = 1; ; ++k) {
                const int x = k * x0;
                const int y = k * y0;
                if (x > limit || y > limit) {
                    break;
                }

                adj[static_cast<std::size_t>(x)].push_back(y);
                adj[static_cast<std::size_t>(y)].push_back(x);
            }
        }
    }

    for (auto& v : adj) {
        std::sort(v.begin(), v.end());
        v.erase(std::unique(v.begin(), v.end()), v.end());
    }

    std::set<int> sums;

    for (int p = 1; p <= limit; ++p) {
        const auto& vp = adj[static_cast<std::size_t>(p)];
        if (vp.size() < 2) {
            continue;
        }

        for (int q : vp) {
            if (q <= p) {
                continue;
            }
            if (p + q >= limit) {
                continue;
            }

            const auto& vq = adj[static_cast<std::size_t>(q)];

            auto it_p = std::upper_bound(vp.begin(), vp.end(), q);
            auto it_q = std::upper_bound(vq.begin(), vq.end(), q);

            while (it_p != vp.end() && it_q != vq.end()) {
                if (*it_p == *it_q) {
                    const int r = *it_p;
                    const int s = p + q + r;
                    if (s <= limit) {
                        sums.insert(s);
                    }
                    ++it_p;
                    ++it_q;
                } else if (*it_p < *it_q) {
                    ++it_p;
                } else {
                    ++it_q;
                }
            }
        }
    }

    return sums;
}

u64 solve(const int limit) {
    const std::set<int> sums = generate_distinct_sums(limit);
    u64 total = 0;
    for (int s : sums) {
        total += static_cast<u64>(s);
    }
    return total;
}

u64 brute_small(const int limit) {
    auto is_square = [](const int x) {
        const int r = static_cast<int>(std::sqrt(static_cast<long double>(x)));
        return r * r == x || (r + 1) * (r + 1) == x;
    };

    std::set<int> sums;
    for (int p = 1; p <= limit; ++p) {
        for (int q = p + 1; p + q <= limit; ++q) {
            if (!is_square(p * p + p * q + q * q)) {
                continue;
            }
            for (int r = q + 1; p + q + r <= limit; ++r) {
                if (is_square(p * p + p * r + r * r) && is_square(q * q + q * r + r * r)) {
                    sums.insert(p + q + r);
                }
            }
        }
    }

    u64 total = 0;
    for (int s : sums) {
        total += static_cast<u64>(s);
    }
    return total;
}

bool run_checkpoints() {
    const std::set<int> sums_1000 = generate_distinct_sums(1000);
    if (sums_1000.count(784) == 0) {
        std::cerr << "Checkpoint failed: 784 not generated for limit 1000" << '\n';
        return false;
    }
    if (solve(1000) != 784ULL) {
        std::cerr << "Checkpoint failed for limit 1000" << '\n';
        return false;
    }
    if (solve(1500) != brute_small(1500)) {
        std::cerr << "Checkpoint failed for brute-force cross-check limit 1500" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << solve(options.limit) << '\n';
    return 0;
}

Python

import math
from typing import Set, List

def generate_distinct_sums(limit: int) -> Set[int]:
    adj = [[] for _ in range(limit + 1)]
    
    m_max = int(math.isqrt(2 * limit)) + 4
    for m in range(2, m_max + 1):
        for n in range(1, m):
            if math.gcd(m, n) != 1:
                continue
            if (m - n) % 3 == 0:
                continue
            
            x0 = m * m - n * n
            y0 = n * (2 * m + n)
            if x0 <= 0 or y0 <= 0:
                continue
            
            if x0 > limit or y0 > limit:
                continue
            
            k = 1
            while True:
                x = k * x0
                y = k * y0
                if x > limit or y > limit:
                    break
                
                adj[x].append(y)
                adj[y].append(x)
                k += 1
    
    for i in range(len(adj)):
        adj[i].sort()
        adj[i] = list(dict.fromkeys(adj[i]))  # Remove duplicates while preserving order
    
    sums: Set[int] = set()
    
    for p in range(1, limit + 1):
        vp = adj[p]
        if len(vp) < 2:
            continue
        
        for q in vp:
            if q <= p:
                continue
            if p + q >= limit:
                continue
            
            vq = adj[q]
            
            # Find first element > q in both lists
            it_p_idx = 0
            while it_p_idx < len(vp) and vp[it_p_idx] <= q:
                it_p_idx += 1
            
            it_q_idx = 0
            while it_q_idx < len(vq) and vq[it_q_idx] <= q:
                it_q_idx += 1
            
            while it_p_idx < len(vp) and it_q_idx < len(vq):
                if vp[it_p_idx] == vq[it_q_idx]:
                    r = vp[it_p_idx]
                    s = p + q + r
                    if s <= limit:
                        sums.add(s)
                    it_p_idx += 1
                    it_q_idx += 1
                elif vp[it_p_idx] < vq[it_q_idx]:
                    it_p_idx += 1
                else:
                    it_q_idx += 1
    
    return sums

def solve(limit: int) -> int:
    sums = generate_distinct_sums(limit)
    return sum(sums)

def brute_small(limit: int) -> int:
    def is_square(x: int) -> bool:
        r = int(math.isqrt(x))
        return r * r == x or (r + 1) * (r + 1) == x
    
    sums: Set[int] = set()
    
    for p in range(1, limit + 1):
        for q in range(p + 1, limit - p + 1):
            if not is_square(p * p + p * q + q * q):
                continue
            for r in range(q + 1, limit - p - q + 1):
                if is_square(p * p + p * r + r * r) and is_square(q * q + q * r + r * r):
                    sums.add(p + q + r)
    
    return sum(sums)

def run_checkpoints() -> bool:
    sums_1000 = generate_distinct_sums(1000)
    if 784 not in sums_1000:
        print("Checkpoint failed: 784 not generated for limit 1000", file=__import__('sys').stderr)
        return False
    if solve(1000) != 784:
        print("Checkpoint failed for limit 1000", file=__import__('sys').stderr)
        return False
    if solve(1500) != brute_small(1500):
        print("Checkpoint failed for brute-force cross-check limit 1500", file=__import__('sys').stderr)
        return False
    return True

def main():
    import sys
    
    limit = 120000
    run_checkpoints_flag = True
    
    i = 1
    while i < len(sys.argv):
        arg = sys.argv[i]
        if arg == "--skip-checkpoints":
            run_checkpoints_flag = False
        elif arg.startswith("--limit="):
            try:
                limit = int(arg[8:])
            except ValueError:
                print(f"Unknown argument: {arg}", file=sys.stderr)
                sys.exit(1)
        else:
            print(f"Unknown argument: {arg}", file=sys.stderr)
            sys.exit(1)
        i += 1
    
    if limit < 3:
        print("Limit must be at least 3", file=sys.stderr)
        sys.exit(1)
    
    if run_checkpoints_flag and not run_checkpoints():
        sys.exit(2)
    
    print(solve(limit))

if __name__ == "__main__":
    main()

Java

import java.util.*;

public class Euler143 {
    static boolean isSquare(long n) {
        long s = (long) Math.sqrt(n);
        return s * s == n;
    }

    public static void main(String[] args) {
        int limit = 120000;
        Map<Integer, Set<Integer>> pairs = new HashMap<>();
        for (int p = 1; p <= limit; p++)
            for (int q = p + 1; p + q <= limit; q++) {
                if (isSquare((long) p * p + (long) p * q + (long) q * q)) {
                    pairs.computeIfAbsent(p, k -> new TreeSet<>()).add(q);
                    pairs.computeIfAbsent(q, k -> new TreeSet<>()).add(p);
                }
            }
        Set<Integer> sums = new TreeSet<>();
        for (int p = 1; p <= limit; p++) {
            if (!pairs.containsKey(p))
                continue;
            Integer[] qs = pairs.get(p).toArray(new Integer[0]);
            for (int i = 0; i < qs.length; i++) {
                int q = qs[i];
                if (p + q >= limit)
                    break;
                for (int j = i + 1; j < qs.length; j++) {
                    int r = qs[j];
                    if (p + q + r > limit)
                        break;
                    Set<Integer> qp = pairs.get(q);
                    if (qp != null && qp.contains(r))
                        sums.add(p + q + r);
                }
            }
        }
        long total = 0;
        for (int s : sums)
            total += s;
        System.out.println(total);
    }
}