Problem 264: Triangle Centres
View on Project EulerProject Euler Problem 264 Solution
EulerSolve provides an optimized solution for Project Euler Problem 264, Triangle Centres, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary We search for triangles with integer vertices \(A,B,C\in\mathbb{Z}^2\) such that all three points lie on the same circle centered at the origin and the centroid is fixed by $$A+B+C=(5,0).$$ For every such triangle whose perimeter does not exceed the given limit, we add its perimeter to the total. The solver prints the final sum rounded to four decimal places. Mathematical Approach 1. Turn the centroid condition into a vector equation Because the centroid is \((5/3,0)\), the vertex sum is the integer vector $$H=(5,0).$$ If we choose one vertex \(C=(c_x,c_y)\), then the remaining two vertices satisfy $$A+B=H-C=(5-c_x,-c_y).$$ The code calls this vector \(S\). Once \(C\) is fixed, the whole problem becomes: find integer pairs \(A,B\) with sum \(S\), equal distance from the origin, and correct perimeter. 2. Split the pair \(A,B\) into sum and difference Write $$A=\frac{S+U}{2},\qquad B=\frac{S-U}{2},\qquad U:=A-B.$$ Then \(A\) and \(B\) have equal radius if and only if $$|A|^2=|B|^2 \iff U\perp S.$$ Moreover, if the common radius is \(R\), then $$|U|^2=4R^2-|S|^2,$$ because $$4|A|^2=|S+U|^2=|S|^2+|U|^2,$$ and the same expression holds for \(B\). This is the central identity used by the code. 3. Integer orthogonal vectors via a gcd basis Since \(S=(s_x,s_y)\) is integral, every integer vector orthogonal to \(S\) lies on a one-dimensional lattice....
Detailed mathematical approach
Problem Summary
We search for triangles with integer vertices \(A,B,C\in\mathbb{Z}^2\) such that all three points lie on the same circle centered at the origin and the centroid is fixed by
$$A+B+C=(5,0).$$
For every such triangle whose perimeter does not exceed the given limit, we add its perimeter to the total. The solver prints the final sum rounded to four decimal places.
Mathematical Approach
1. Turn the centroid condition into a vector equation
Because the centroid is \((5/3,0)\), the vertex sum is the integer vector
$$H=(5,0).$$
If we choose one vertex \(C=(c_x,c_y)\), then the remaining two vertices satisfy
$$A+B=H-C=(5-c_x,-c_y).$$
The code calls this vector \(S\). Once \(C\) is fixed, the whole problem becomes: find integer pairs \(A,B\) with sum \(S\), equal distance from the origin, and correct perimeter.
2. Split the pair \(A,B\) into sum and difference
Write
$$A=\frac{S+U}{2},\qquad B=\frac{S-U}{2},\qquad U:=A-B.$$
Then \(A\) and \(B\) have equal radius if and only if
$$|A|^2=|B|^2 \iff U\perp S.$$
Moreover, if the common radius is \(R\), then
$$|U|^2=4R^2-|S|^2,$$
because
$$4|A|^2=|S+U|^2=|S|^2+|U|^2,$$
and the same expression holds for \(B\). This is the central identity used by the code.
3. Integer orthogonal vectors via a gcd basis
Since \(S=(s_x,s_y)\) is integral, every integer vector orthogonal to \(S\) lies on a one-dimensional lattice. Let
$$g=\gcd(|s_x|,|s_y|),\qquad P=\left(-\frac{s_y}{g},\frac{s_x}{g}\right).$$
Then \(P\) is the primitive integer normal vector to \(S\), and every integer orthogonal vector is
$$U=tP,\qquad t\in\mathbb{Z}.$$
Substituting this into the radius identity gives
$$t^2=\frac{(4R^2-|S|^2)g^2}{|S|^2}.$$
The solver checks exactly this quantity: it must be an integer square before \(A\) and \(B\) are even constructed.
4. Parity, degeneracy, and canonical uniqueness
Once \(U\) is known, the coordinates of \(A\) and \(B\) are obtained from \((S\pm U)/2\). Therefore each coordinate must have the correct parity; otherwise the pair is discarded immediately.
The code then rejects degenerate triangles by checking that the doubled area, computed by a cross product, is nonzero.
To avoid counting the same triangle multiple times, the code keeps only a canonical representative: the chosen vertex \(C\) must be lexicographically no larger than \(A\) and \(B\). This removes the six permutations of the same unlabeled triangle, while the fixed sign choice for the primitive normal vector removes the \(A\leftrightarrow B\) swap.
5. Why the search for \(C\) is bounded
The code does not scan all possible radii blindly. It first derives an upper bound on the circumradius from the perimeter cap.
Using
$$AB^2=4R^2-|A+B|^2=4R^2-|H-C|^2,$$
and the triangle inequality \(|H-C|\le |H|+|C|=5+R\), we get
$$AB^2\ge 4R^2-(R+5)^2=3R^2-10R-25.$$
The same lower bound applies symmetrically to the three sides, so
$$P\ge 3\sqrt{3R^2-10R-25}.$$
Solving this inequality for \(R\) gives the radius cutoff used in the program. After that, the solver only scans integer points \(C=(c_x,c_y)\) with \(|C|\le R_{\max}\), and because \(C\) is the lexicographically smallest vertex, its \(x\)-coordinate must satisfy \(c_x\le 1\). That is why the outer scan stops at \(x=1\).
6. Threading and numeric robustness
The implementation splits the \(c_x\)-range into chunks and distributes them across threads. Each thread keeps a local perimeter sum and triangle count, and the partial sums are merged at the end.
Geometric tests that need exactness use integer arithmetic and the integer square root helper. The final perimeter is computed with long double, then rounded to four decimals. The checkpoints in the code verify both a sample perimeter value and that the multi-threaded and single-threaded runs agree.
The checkpoints are concrete: the code requires round4(solve(50, 1)) = 291.0089, and it also checks that solve(300, 1) matches a multi-threaded run. One checkpoint validates the geometric pipeline numerically, and the other validates that the thread split does not change the result.
7. Worked structural example
Take the valid triangle
$$C=(-4,-3),\qquad A=(4,3),\qquad B=(5,0).$$
All three vertices satisfy
$$|A|^2=|B|^2=|C|^2=25,$$
and also
$$A+B+C=(4+5-4,\ 3+0-3)=(5,0).$$
For this choice of \(C\), we get
$$S=H-C=(9,3),\qquad g=\gcd(9,3)=3,\qquad P=(-1,3).$$
The radius identity gives
$$t^2=\frac{(4\cdot 25-90)\cdot 3^2}{90}=1,$$
so \(t=1\) and therefore
$$U=tP=(-1,3).$$
Reconstructing the pair yields
$$A=\frac{S+U}{2}=(4,3),\qquad B=\frac{S-U}{2}=(5,0).$$
This is exactly what the code does for every admissible \(C\): derive \(S\), test whether the induced \(t^2\) is a perfect square, then rebuild \(A\) and \(B\) algebraically instead of enumerating all vertex triples.
How the Code Works
The program accepts --perimeter, --threads, and --skip-checkpoints. The function solve first converts the perimeter limit into a radius bound \(R_{\max}\). It then scans integer points \(C=(c_x,c_y)\) inside that circle. For each \(C\), it forms
$$S=(5-c_x,-c_y),\qquad r^2=|C|^2.$$
After that it computes the primitive orthogonal vector \(P\), checks that the derived \(t^2\) is a perfect square, verifies parity, reconstructs \(A\) and \(B\), rejects degenerate or non-canonical triangles, and finally filters by the exact perimeter bound.
The helper isqrt_i64 is used both for radius estimates and square tests. The helper lex_leq implements the canonical ordering. The checkpoint routine compares the single-threaded and multi-threaded results and also checks the stored sample sum for perimeter \(50\).
Complexity Analysis
The dominant work is the scan over candidate vertices \(C\) in the bounded disk. For each candidate, the solver performs only constant-time arithmetic, gcd, square-root, and geometry checks. Because the scan is restricted to \(c_x\le 1\) and \(|C|\le R_{\max}\), the search space is much smaller than all triples of vertices.
Time grows roughly quadratically with the radius bound, and the multi-threaded split reduces wall-clock time without changing the mathematics. Memory usage stays small because the solver keeps only a few counters and partial sums.
Footnotes and References
- Problem page: https://projecteuler.net/problem=264
- Primitive lattice vectors and Gaussian integers: Wikipedia - Gaussian integer
- Sum/difference decomposition: Wikipedia - Parallelogram law
- Integer square root and exact arithmetic: Wikipedia - Integer square root
Problem 264 source code
C++
#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <string>
#include <thread>
#include <vector>
namespace {
using i64 = long long;
using u64 = std::uint64_t;
struct Options {
int perimeter_limit = 100'000;
int threads = static_cast<int>(std::thread::hardware_concurrency());
bool run_checkpoints = true;
};
struct PointI {
i64 x = 0;
i64 y = 0;
};
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) {
if (options.threads <= 0) {
options.threads = 1;
}
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, "--perimeter=", options.perimeter_limit) ||
parse_int_after_prefix(arg, "--threads=", options.threads)) {
continue;
}
std::cerr << "Unknown argument: " << arg << '\n';
return false;
}
return options.perimeter_limit >= 1 && options.threads >= 1;
}
bool lex_leq(const PointI& a, const PointI& b) {
if (a.x != b.x) {
return a.x < b.x;
}
return a.y <= b.y;
}
bool point_equal(const PointI& a, const PointI& b) {
return a.x == b.x && a.y == b.y;
}
i64 isqrt_i64(i64 x) {
i64 r = static_cast<i64>(std::sqrt(static_cast<long double>(x)));
while ((r + 1) * (r + 1) <= x) {
++r;
}
while (r * r > x) {
--r;
}
return r;
}
long double dist(const PointI& a, const PointI& b) {
const long double dx = static_cast<long double>(a.x - b.x);
const long double dy = static_cast<long double>(a.y - b.y);
return std::sqrt(dx * dx + dy * dy);
}
long double solve(const int perimeter_limit, int threads) {
// From AB^2 = 4R^2 - |A+B|^2 with A+B = H-C and |H|=5,
// each side has a lower bound sqrt(3R^2 - 10R - 25).
// So p >= 3*sqrt(3R^2 - 10R - 25).
const long double p = static_cast<long double>(perimeter_limit);
const long double r_bound =
(10.0L + std::sqrt(100.0L + 12.0L * (25.0L + (p * p) / 9.0L))) / 6.0L;
const int r_max = static_cast<int>(std::floor(r_bound));
threads = std::max(1, std::min(threads, r_max + 2));
std::atomic<int> next_x{-r_max};
constexpr int kChunk = 32;
std::vector<long double> partial_sum(static_cast<std::size_t>(threads), 0.0L);
std::vector<u64> partial_count(static_cast<std::size_t>(threads), 0);
std::vector<std::thread> pool;
pool.reserve(static_cast<std::size_t>(threads));
for (int t = 0; t < threads; ++t) {
pool.emplace_back([&, t]() {
long double local_sum = 0.0L;
u64 local_count = 0;
while (true) {
const int x_start = next_x.fetch_add(kChunk, std::memory_order_relaxed);
if (x_start > 1) {
break;
}
const int x_end = std::min(1, x_start + kChunk - 1);
for (int cx = x_start; cx <= x_end; ++cx) {
const i64 cx2 = static_cast<i64>(cx) * static_cast<i64>(cx);
const i64 cy_max = isqrt_i64(static_cast<i64>(r_max) * r_max - cx2);
for (i64 cy = -cy_max; cy <= cy_max; ++cy) {
const i64 r2 = cx2 + cy * cy;
if (r2 == 0) {
continue;
}
const i64 sx = 5 - static_cast<i64>(cx);
const i64 sy = -cy;
const i64 den = sx * sx + sy * sy;
if (den == 0) {
continue;
}
const i64 d = 4 * r2 - den;
if (d <= 0) {
continue;
}
const i64 g = std::gcd(std::llabs(sx), std::llabs(sy));
const i64 num = d * g * g;
if (num % den != 0) {
continue;
}
const i64 t2 = num / den;
const i64 tt = isqrt_i64(t2);
if (tt * tt != t2) {
continue;
}
const i64 px = -sy / g;
const i64 py = sx / g;
const i64 ux = tt * px;
const i64 uy = tt * py;
if (((sx + ux) & 1LL) != 0LL || ((sy + uy) & 1LL) != 0LL) {
continue;
}
const PointI a{(sx + ux) / 2, (sy + uy) / 2};
const PointI b{(sx - ux) / 2, (sy - uy) / 2};
const PointI c{cx, cy};
if (point_equal(a, b) || point_equal(a, c) || point_equal(b, c)) {
continue;
}
if (!lex_leq(c, a) || !lex_leq(c, b)) {
continue;
}
if (a.x * a.x + a.y * a.y != r2 || b.x * b.x + b.y * b.y != r2) {
continue;
}
if (a.x + b.x + c.x != 5 || a.y + b.y + c.y != 0) {
continue;
}
const i64 area2 = std::llabs((b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x));
if (area2 == 0) {
continue;
}
const long double perimeter = dist(a, b) + dist(a, c) + dist(b, c);
if (perimeter > static_cast<long double>(perimeter_limit) + 1e-12L) {
continue;
}
++local_count;
local_sum += perimeter;
}
}
}
partial_sum[static_cast<std::size_t>(t)] = local_sum;
partial_count[static_cast<std::size_t>(t)] = local_count;
});
}
for (auto& th : pool) {
th.join();
}
long double total = 0.0L;
u64 count = 0;
for (int t = 0; t < threads; ++t) {
total += partial_sum[static_cast<std::size_t>(t)];
count += partial_count[static_cast<std::size_t>(t)];
}
(void)count;
return total;
}
long double round4(long double x) {
return std::round(x * 10000.0L) / 10000.0L;
}
bool run_checkpoints() {
const long double sum50 = round4(solve(50, 1));
if (std::fabsl(sum50 - 291.0089L) > 1e-9L) {
std::cerr << "Checkpoint failed for perimeter=50: got " << std::setprecision(10)
<< static_cast<double>(sum50) << '\n';
return false;
}
unsigned hw = std::thread::hardware_concurrency();
if (hw == 0) {
hw = 2;
}
const int multi_threads = static_cast<int>(std::min<unsigned>(hw, 8));
const long double s1 = round4(solve(300, 1));
const long double sm = round4(solve(300, multi_threads));
if (std::fabsl(s1 - sm) > 1e-9L) {
std::cerr << "Thread consistency failed for perimeter=300" << '\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;
}
const long double answer = round4(solve(options.perimeter_limit, options.threads));
std::cout << std::fixed << std::setprecision(4) << static_cast<double>(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 ""
answer_candidates = []
equal_candidates = []
for line in lines:
m1 = ANSWER_RE.search(line)
if m1:
answer_candidates.append(m1.group(1).strip())
m2 = EQUAL_RE.search(line)
if m2:
equal_candidates.append(m2.group(1).strip())
if answer_candidates:
return answer_candidates[-1]
if equal_candidates:
return equal_candidates[-1]
return lines[-1]
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 = subprocess.check_output([str(binary)], text=True)
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.util.*;
import java.util.concurrent.*;
public class Euler264 {
static long isqrt(long x) {
if (x < 0)
return 0;
long r = (long) Math.sqrt(x);
while ((r + 1) * (r + 1) <= x)
r++;
while (r * r > x)
r--;
return r;
}
static double dist(long ax, long ay, long bx, long by) {
double dx = ax - bx;
double dy = ay - by;
return Math.sqrt(dx * dx + dy * dy);
}
static long gcd(long a, long b) {
while (b != 0) {
long t = b;
b = a % b;
a = t;
}
return Math.abs(a);
}
public static String solve() {
int perimeterLimit = 100000;
int threads = Math.max(1, Runtime.getRuntime().availableProcessors());
double p = perimeterLimit;
double rBound = (10.0 + Math.sqrt(100.0 + 12.0 * (25.0 + (p * p) / 9.0))) / 6.0;
int rMax = (int) Math.floor(rBound);
threads = Math.max(1, Math.min(threads, rMax + 2));
int chunk = 32;
List<int[]> tasks = new ArrayList<>();
int startX = -rMax;
while (startX <= 1) {
int endX = Math.min(1, startX + chunk - 1);
tasks.add(new int[] { startX, endX });
startX += chunk;
}
double totalSum = 0.0;
ExecutorService executor = Executors.newFixedThreadPool(threads);
List<Future<Double>> futures = new ArrayList<>();
for (int[] task : tasks) {
final int xStart = task[0];
final int xEnd = task[1];
futures.add(executor.submit(() -> {
double localSum = 0.0;
for (int cx = xStart; cx <= xEnd; ++cx) {
long cx2 = (long) cx * cx;
long cyMax = isqrt((long) rMax * rMax - cx2);
for (long cy = -cyMax; cy <= cyMax; ++cy) {
long r2 = cx2 + cy * cy;
if (r2 == 0)
continue;
long sx = 5 - cx;
long sy = -cy;
long den = sx * sx + sy * sy;
if (den == 0)
continue;
long d = 4 * r2 - den;
if (d <= 0)
continue;
long g = gcd(Math.abs(sx), Math.abs(sy));
long num = d * g * g;
if (num % den != 0)
continue;
long t2 = num / den;
long tt = isqrt(t2);
if (tt * tt != t2)
continue;
long px = -sy / g;
long py = sx / g;
long ux = tt * px;
long uy = tt * py;
if (((sx + ux) & 1) != 0 || ((sy + uy) & 1) != 0)
continue;
long ax = (sx + ux) / 2;
long ay = (sy + uy) / 2;
long bx = (sx - ux) / 2;
long by = (sy - uy) / 2;
if ((ax == bx && ay == by) || (ax == cx && ay == cy) || (bx == cx && by == cy))
continue;
if (!(cx < ax || (cx == ax && cy <= ay)))
continue;
if (!(cx < bx || (cx == bx && cy <= by)))
continue;
if (ax * ax + ay * ay != r2 || bx * bx + by * by != r2)
continue;
if (ax + bx + cx != 5 || ay + by + cy != 0)
continue;
long area2 = Math.abs((bx - ax) * (cy - ay) - (by - ay) * (cx - ax));
if (area2 == 0)
continue;
double perimeter = dist(ax, ay, bx, by) + dist(ax, ay, cx, cy) + dist(bx, by, cx, cy);
if (perimeter > perimeterLimit + 1e-12)
continue;
localSum += perimeter;
}
}
return localSum;
}));
}
for (Future<Double> f : futures) {
try {
totalSum += f.get();
} catch (Exception e) {
}
}
executor.shutdown();
double roundTotal = Math.round(totalSum * 10000.0) / 10000.0;
return String.format(Locale.US, "%.4f", roundTotal);
}
public static void main(String[] args) {
System.out.println(solve());
}
}