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 $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.

The basic arithmetic object: a 120-degree integer pair

Call two positive integers $x$ and $y$ compatible if

$$x^2+xy+y^2=z^2$$

for some integer $z$. Geometrically, $x$ and $y$ can then appear as two Torricelli distances meeting at $120^\circ$, and $z$ is the opposite triangle side. The full problem is therefore not about arbitrary triples first; it is about finding three distances that are pairwise compatible.

That observation turns the search into a graph problem. Take the integers $1,2,\dots,N$ with $N=120{,}000$ as vertices, and join $x$ and $y$ by an edge whenever they are compatible. Then every valid $(p,q,r)$ is exactly a 3-clique in this graph, because all three pairwise square conditions must hold simultaneously.

Generating all compatible pairs efficiently

The optimized implementations do not test every pair $(x,y)$ from scratch. They use the standard parametrization of primitive solutions of

$$x^2+xy+y^2=z^2,$$

namely, up to swapping $x$ and $y$,

$$x_0=m^2-n^2,\qquad y_0=n(2m+n),\qquad z_0=m^2+mn+n^2,$$

with

$$m>n>0,\qquad \gcd(m,n)=1,\qquad m-n\not\equiv 0 \pmod 3.$$

Every non-primitive compatible pair is then just a multiple $(kx_0,ky_0)$. This is why the fast versions loop over $(m,n)$, generate primitive seeds, and then scale them by $k$ until the limit is reached. The side length $z$ never needs to be stored, because the existence of an edge is all that matters later.

Why triangle enumeration becomes common-neighbor intersection

Suppose $p<q<r$. In graph language, $(p,q,r)$ is valid exactly when $q$ and $r$ are both neighbors of $p$, and $r$ is also a neighbor of $q$. So after building sorted adjacency lists, one can fix $p$, choose $q>p$ among the neighbors of $p$, and then look for common neighbors of $p$ and $q$ that are larger than $q$. Intersecting the two sorted tails finds precisely those $r$.

This ordered search has two important invariants. First, $p<q<r$ prevents the same triple from being generated in six different permutations. Second, the problem asks for distinct sums $p+q+r$, not distinct triples, so the final sums are inserted into a set before being added together. Two different triples may contribute the same perimeter, and that perimeter must be counted once.

Worked example: the perimeter 784

A concrete solution is

$$ (p,q,r)=(195,264,325). $$

Indeed,

$$195^2+195\cdot264+264^2=399^2,$$

$$195^2+195\cdot325+325^2=455^2,$$

$$264^2+264\cdot325+325^2=511^2.$$

So the three sides of the triangle are $399$, $455$, and $511$, all integral, and the corresponding Torricelli-distance sum is

$$195+264+325=784.$$

This is the small example used by the implementations as a sanity check before attacking the full limit.

How the Code Works

Building the compatibility graph

The C++, Python, and Java implementations all work with the same graph of compatible distance pairs. The optimized C++ and Python versions build that graph from the parametrization above, adding each generated pair to both adjacency lists because compatibility is symmetric. After generation, each list is sorted and deduplicated so that later intersections are linear scans over ordered data.

Closing 3-cliques into Torricelli triples

Once the graph exists, the next phase is to enumerate every ordered triple $p<q<r$ that forms a clique. For a fixed $p$, the implementation scans its neighbors $q$ with $q>p$. Any admissible $r$ must be a common neighbor of $p$ and $q$ and must also satisfy $r>q$. The sorted-list intersection enforces all of that at once, and whenever $p+q+r \le 120{,}000$, the sum is recorded.

Optimized versus direct construction

The Java implementation illustrates the same mathematics in a more direct style. Instead of generating compatible pairs from the parametrization, it checks every $p<q$ with $p+q \le N$, tests whether $p^2+pq+q^2$ is a square, stores the successful pairs, and then searches for common neighbors to complete triangles. That approach is conceptually simple but spends much more time on failed pair tests. The C++ and Python implementations avoid most of that work by generating only genuine compatible pairs from the start.

Complexity Analysis

Let $E$ be the number of compatible pairs up to the limit $N$. The optimized construction spends $O(E)$ time generating scaled pairs from primitive seeds, plus the cost of sorting the adjacency lists, plus the cost of all ordered list intersections used to enumerate cliques. A conservative worst-case bound is still quadratic in $N$, but the arithmetic graph is sparse enough that the full limit is practical.

The direct construction used in the Java version performs $O(N^2)$ pair tests just to build the graph, and only afterward begins the common-neighbor search. All versions use $O(E)$ memory for the adjacency structure, together with a set of distinct perimeter sums. The main speedup comes from replacing a naive search over triples of distances with a search over edges and their common neighbors.

Footnotes and References

Problem page: Project Euler 143
Torricelli or Fermat point: Wikipedia - Fermat point
Cosine law: Wikipedia - Law of cosines
The quadratic form $x^2+xy+y^2$ as an Eisenstein norm: Wikipedia - Eisenstein integer

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);
    }
}