Problem 264: Triangle Centres

Project 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 look at triangles $ABC$ whose three vertices are lattice points, that is, points with integer coordinates, whose circumcentre, the centre of the circle through $A$, $B$ and $C$, is the origin $O$, and whose orthocentre, the common point of the three altitudes, is $H=(5,0)$. The statement lists the nine such triangles with perimeter at most $50$; the sum of their perimeters, rounded to four decimal places, is $291.0089$. We must find all such triangles with perimeter at most $10^5$ and give the sum of their perimeters, rounded to four decimal places. ^[1]

The whole solution in one sentence: when the circumcentre is the origin, the orthocentre is simply $A+B+C$, so a single vertex $C$ other than $H$ already fixes the opposite side as the chord of the circle with midpoint $(H-C)/2$, and the program scans the lattice points $C$ of a disc whose radius follows from the perimeter limit and keeps those for which this chord ends at lattice points.

Mathematical Approach

We treat each point $P=(p_x,p_y)$ as the vector from the origin to it, with length $|P|=\sqrt{p_x^2+p_y^2}$. The dot product of $P$ and $Q=(q_x,q_y)$ is $P\cdot Q=p_xq_x+p_yq_y$. Expanding gives $|P-Q|^2=|P|^2+|Q|^2-2\,P\cdot Q$, so by Pythagoras' theorem and its converse $P\cdot Q=0$ exactly when $P$ and $Q$ are perpendicular, or one of them is zero; for example $(1,0)\cdot(0,1)=0$. $R$ denotes the circumradius, so $|A|=|B|=|C|=R$, and $R^2$ is an integer. The triangle $A=(1,8)$, $B=(8,-1)$, $C=(-4,-7)$ from the statement's list, with $R^2=65$, serves as the running example.

1. The orthocentre is the sum of the vertices

Put $H'=A+B+C$. Then

$$\begin{aligned}(H'-A)\cdot(B-C)&=(B+C)\cdot(B-C)\\&=|B|^2-|C|^2=0,\end{aligned}$$

so $H'$ lies on the line through $A$ perpendicular to $BC$, the altitude from $A$. In the same way $H'$ lies on the altitudes from $B$ and $C$. Two altitudes of a triangle meet in exactly one point, so $H'$ is the orthocentre; this is Sylvester's relation $\vec{OH}=\vec{OA}+\vec{OB}+\vec{OC}$. ^[2] The problem therefore asks for the lattice points with

$$|A|=|B|=|C|=R,\qquad A+B+C=H.$$

The running example passes: $(1,8)+(8,-1)+(-4,-7)=(5,0)$ and $1+64=64+1=16+49=65$. The reflection $(x,y)\mapsto(x,-y)$ fixes $O$ and $H$, so the solutions come in mirror pairs; of the nine listed triangles only the one with the vertices $(0,5)$, $(5,0)$, $(0,-5)$ is its own mirror image.

2. Small cases: the five right triangles

The first five triangles of the list share the vertex $H$ itself. If $C=H$, then $A+B=0$, so $AB$ is a diameter, and by Thales's theorem the angle at $C$ is a right angle. ^[3] As $H$ then lies on the circle, $R=5$. The circle $x^2+y^2=25$ contains exactly twelve lattice points, $(\pm5,0)$, $(0,\pm5)$, $(\pm3,\pm4)$ and $(\pm4,\pm3)$. Without $(\pm5,0)$, for which $A$ or $B$ would coincide with $C$, five antipodal pairs remain, and they give five right triangles with hypotenuse $10$ and the perimeters $10+4\sqrt{10}\approx22.6491$ (twice), $10+6\sqrt5\approx23.4164$ (twice) and $10+10\sqrt2\approx24.1421$. For instance $A=(4,3)$ and $B=(-4,-3)$ give the legs $|AH|=\sqrt{10}$ and $|BH|=\sqrt{90}=3\sqrt{10}$.

3. A vertex other than $H$ fixes the triangle

Fix a lattice point $C$ with $|C|=R$ and put

$$S=H-C=A+B,\qquad U=A-B.$$

Then $|A|^2-|B|^2=(A+B)\cdot(A-B)$${}=S\cdot U$, and the parallelogram law gives $2|A|^2+2|B|^2=|S|^2+|U|^2$. ^[4] Hence $A$ and $B$ lie on the circle of radius $R$ exactly when

$$S\cdot U=0,\qquad |U|^2=4R^2-|S|^2,$$

and then $A=\frac12(S+U)$ and $B=\frac12(S-U)$. If $C\ne H$, the midpoint $\frac12S$ of $AB$ is not the centre, and geometrically $AB$ is the chord perpendicular to the radius through that midpoint. So $C$ determines $U$ up to its sign, and changing the sign only swaps $A$ and $B$: every point $C\ne H$ is a vertex of at most one valid triangle. The exception $C=H$ gives $S=0$, and then $AB$ can be any of the five diameters of Section 2. The scan of Section 7 only uses vertices with $c_x\le1$, so it never starts from $H$, and it finds those five triangles from their other vertices. Conversely, every such $U$ with $4R^2-|S|^2\gt0$ gives two different points $A$ and $B$ of the circle with $A+B+C=H$. If both also differ from $C$, the three points are not collinear, because a line meets a circle in at most two points, and Section 1 makes $H$ the orthocentre.

4. Lattice vectors perpendicular to $S$

Write $S=(s_x,s_y)\ne(0,0)$, let $g=\gcd(s_x,s_y)$ and

$$P=\Bigl(-\frac{s_y}{g},\ \frac{s_x}{g}\Bigr).$$

The integer vectors $U=(u_x,u_y)$ with $S\cdot U=0$ are exactly the multiples $tP$ with $t\in\mathbb Z$. If $s_x=0$, then $s_y\ne0$, so $S\cdot U=0$ forces $u_y=0$, while $P=(\mp1,0)$; again $U$ is a multiple of $P$. Otherwise, dividing $s_xu_x=-s_yu_y$ by $g$ leaves the coprime coefficients $s_x/g$ and $s_y/g$, so by Euclid's lemma ^[5] $s_x/g$ divides $u_y$; writing $u_y=t\,s_x/g$ and cancelling the nonzero $s_x$ gives $u_x=-t\,s_y/g$. Since $|P|^2=|S|^2/g^2$, the length condition of Section 3 becomes

$$t^2=\frac{(4R^2-|S|^2)\,g^2}{|S|^2}.$$

So $C$ is a vertex of a valid triangle exactly when $4R^2-|S|^2\gt0$, this fraction is an integer, the integer is a perfect square $t^2$, and the points $A$ and $B$ built from $U=tP$ are lattice points different from $C$. The next section shows that the halving in $A=\frac12(S+U)$ never fails.

5. The halving always works

The vector $W=S+U$ satisfies $|W|^2=|S|^2+|U|^2=4R^2$, because $S$ and $U$ are perpendicular. The square of an even number is divisible by $4$, and the square of an odd number leaves the remainder $1$, so a sum of two squares can only be divisible by $4$ when both numbers are even. Hence $W$ has even coordinates, and so has $S-U=W-2U$: whenever $t^2$ is a perfect square, $A$ and $B$ are lattice points.

6. Worked example

For $C=(-4,-7)$ we get $S=(9,7)$, $|S|^2=130$, $g=1$ and $P=(-7,9)$. Here $4R^2-|S|^2=260-130=130$, so $t^2=130\cdot1/130=1$, $t=1$ and $U=(-7,9)$. Then $S+U=(2,16)$ and $S-U=(16,-2)$, so $A=(1,8)$ and $B=(8,-1)$: the running example. Two failures show the tests at work. For $C=(-7,0)$, with $R^2=49$, we get $S=(12,0)$, $g=12$, $P=(0,1)$ and $t^2=52\cdot144/144=52$, which is not a square: the chord runs from $(6,-\sqrt{13})$ to $(6,\sqrt{13})$. For $C=(-6,-1)$, with $R^2=37$, we get $S=(11,1)$, $g=1$ and the fraction $26/122=13/61$, which is not even an integer.

For the limit $50$ of the statement the scan of Section 9 visits the $220$ lattice points $C\ne O$ with $c_x\le1$ and $|C|\le11$. Of these, $184$ satisfy $4R^2\gt|S|^2$, only $38$ of those make the fraction an integer, and $12$ make it a perfect square; these $12$ points produce the nine listed triangles, three of them twice.

7. Counting every triangle once

The lexicographic order compares points first by $x$ and then by $y$: $P$ comes before $Q$ when $p_x\lt q_x$, or when $p_x=q_x$ and $p_y\lt q_y$. ^[6] It is a total order of the plane, so every triangle has exactly one smallest vertex, and the program keeps a triangle only when the scanned point $C$ is its smallest vertex. The smallest vertex has the smallest $x$-coordinate, so $3c_x\le a_x+b_x+c_x=5$, that is $c_x\le1$, and only these columns are scanned. For the smallest vertex, Section 3 leaves exactly one pair $\{A,B\}$, and the program takes only $U=tP$ with $t\ge1$, one of the two signs. Hence every triangle is counted exactly once. In the example the vertex $(1,8)$ also lies in a scanned column and produces the same triangle, but it is rejected, because $(-4,-7)$ comes first.

8. The side lengths and a bound for the radius

Section 3 gives $|AB|^2=|U|^2=4R^2-|S|^2$, and $|S|^2=|H-C|^2=25-10c_x+R^2$, so

$$|AB|^2=3R^2-25+10\,c_x.$$

Each side is therefore fixed by $R$ and the $x$-coordinate of the opposite vertex. In the example $|AB|^2=195-25-40=130$, $|BC|^2=170+10=180$ and $|CA|^2=170+80=250$, so the perimeter is $\sqrt{130}+\sqrt{180}+\sqrt{250}$${}\approx40.6296$. Adding the three formulas and using $a_x+b_x+c_x=5$ gives $|AB|^2+|BC|^2+|CA|^2=9R^2-25$, which is Euler's formula $OH^2=9R^2-(a^2+b^2+c^2)$ for the side lengths $a,b,c$ and $OH=5$; ^[2] indeed $130+180+250=560=9\cdot65-25$.

Every vertex has an $x$-coordinate of at least $-R$, so every squared side is at least $3R^2-10R-25$. The perimeter $p$ is at least three times the shortest side, so $p^2\ge9\,(3R^2-10R-25)$. If $p\le L$, solving this quadratic inequality for $R$ gives

$$R\le R_L=\frac{10+\sqrt{100+12\,(25+L^2/9)}}{6}.$$

For $L=50$ this is $R_L\approx11.85$, and for $L=10^5$ it is $R_L\approx19246.68$.

9. Why the radius may be rounded down

The program scans the points with $|C|\le r_{\max}=\lfloor R_L\rfloor$. This rounding needs a reason: $R^2$ is an integer but not always a perfect square, so a triangle with $r_{\max}\lt R\le R_L$, for instance with $R^2=r_{\max}^2+1$, would be missed, and the bound of Section 8 does not exclude it, since at $R=19246$ it only guarantees $p\gt99996$. A sharper estimate does. For fixed $R\ge5$, the function $f(x)=\sqrt{3R^2-25+10x}$ is concave, so on the interval $[-R,R]$ its graph lies above the chord joining its end points. ^[7] For $x$ in this interval put $\lambda=(x+R)/(2R)$, so that $0\le\lambda\le1$ and $x=(1-\lambda)(-R)+\lambda R$; concavity then says

$$f(x)\ge(1-\lambda)\,f(-R)+\lambda\,f(R).$$

By Section 8 the perimeter is $f(a_x)+f(b_x)+f(c_x)$. For the three vertices the weights $\lambda$ add up to $(3R+5)/(2R)$, because $a_x+b_x+c_x=5$, so adding the three inequalities gives

$$\begin{aligned}p&\ge\tfrac32\bigl(f(R)+f(-R)\bigr)\\&\quad+\tfrac{5}{2R}\bigl(f(R)-f(-R)\bigr).\end{aligned}$$

The second term is positive, so $p\gt\frac32\bigl(f(R)+f(-R)\bigr)$. This lower bound grows with $R$, and at $R=19246$ it is already about $100005.15$. Hence every triangle with $p\le10^5$ has $R\lt19246$ and lies inside the scanned disc. Evaluating the bound at $r_{\max}$ for every integer limit $23\le L\le10^5$ shows in the same way that rounding down is safe for all of them, and for $L\le22$ there is no triangle at all, since the smallest perimeter is $10+4\sqrt{10}\approx22.65$. In the other direction, $(a-b)^2+(b-c)^2+(c-a)^2\ge0$ gives $p^2\le3(a^2+b^2+c^2)$, so $p\le\sqrt{3(9R^2-25)}\lt3\sqrt3R$. Large solutions are thus nearly equilateral, with a perimeter just below $3\sqrt3R$, while the bound of Section 8 is about $3\sqrt3R-5\sqrt3$; this loss of $5/3$ in the radius is more than the rounding takes back.

How the Code Works

The C++ program ^[8] reads the options --perimeter= (default $100000$) and --skip-checkpoints, runs run_checkpoints unless told to skip it, rounds the result of solve with round4 and prints it with four decimals. solve computes $R_L$ of Section 8 in long double, sets r_max to its floor and walks through the columns $c_x=-r_{\max},\dots,1$ of Section 7 in blocks of $32$. In each column, isqrt_i64 gives the largest $|c_y|$ with $c_x^2+c_y^2\le r_{\max}^2$, and every point $C\ne O$ of the column is tested in this order:

den $=|S|^2$ for $S=(5-c_x,-c_y)$, and d $=4R^2-|S|^2$, which must be positive (Section 3); the test den == 0 can never fire, because $s_x=5-c_x\ge4$;
g $=\gcd(|s_x|,|s_y|)$, computed with std::gcd, and num $=d\,g^2$, which must be divisible by den, and the quotient t2, which must be a perfect square tt*tt (Section 4);
$U=t\,P$, the parity of $S+U$ (Section 5), $A=\frac12(S+U)$, $B=\frac12(S-U)$ and three distinct vertices;
lex_leq(c, a) and lex_leq(c, b), the smallest-vertex rule of Section 7;
re-checks of $|A|^2=|B|^2=R^2$, of $A+B+C=H$ and of a nonzero doubled area, the cross product of $B-A$ and $C-A$;
the perimeter, the sum of three dist values in long double, which is added to the total if it does not exceed $L+10^{-12}$.

These extra tests never reject a point of the scanned columns: the parity test by Section 5, the radius and sum re-checks by Section 3, the distinctness test because $A=C$ would force $B=H-2C$, and $|H-2C|^2=R^2$ means $20c_x=3R^2+25$, so $c_x\ge2$ (likewise for $B$), and the area test because three distinct points of a circle are never collinear. They are safety checks. isqrt_i64 takes the long double square root and corrects it with two loops until $r^2\le x\lt(r+1)^2$, so it returns exactly $\lfloor\sqrt x\rfloor$. ^[9] All integers fit easily into 64 bits: $R^2\le r_{\max}^2\lt3.71\cdot10^8$, and since $g^2\le|S|^2$ and a product of two numbers with the sum $4R^2$ is at most $4R^4$, we get $d\,g^2\le(4R^2-|S|^2)\,|S|^2$${}\le4R^4\lt5.5\cdot10^{17}$, far below $2^{63}\approx9.2\cdot10^{18}$. run_checkpoints requires that the limit $50$ gives $291.0089$ after round4; it also evaluates the limit $300$ twice, once with a single running sum and once with the blocks spread over up to eight partial sums, and requires both rounded results to agree.

The Java program performs the same scan with the same tests, using long integers, its own gcd and isqrt, and double perimeters; the sums of the blocks of $32$ columns are added at the end, and the total is rounded with Math.round and printed with String.format. An exception while a block sum is collected is caught and ignored, so a failing block would silently be missing from the total; normal runs raise none. Its limit is fixed at $100000$, and it runs no checkpoints. The Python file is a bridge: it compiles the C++ source Euler264.cpp from the neighbouring folder solutionsCpp with clang++ or g++ whenever the binary is missing or older than the source, runs it without options, so that the checkpoints run first, and returns the printed line. With GCC on x86, long double uses the 80-bit extended format, but on 64-bit ARM under macOS and Windows it has the same double-precision format as double, ^[10] so the accuracy analysis below assumes double precision.

Complexity and Verification

The scan visits about $\frac\pi2r_{\max}^2$ points, exactly $581\,893\,845$ for $L=10^5$, and spends a bounded number of arithmetic operations and at most one gcd of $O(\log r_{\max})$ steps on each. Since $R_L$ grows linearly with $L$, the time grows like $L^2\log L$. The C++ program needs only a handful of variables; the Java program also keeps two lists with one entry per block of $32$ columns, $602$ entries each for $L=10^5$. Almost every point fails at the divisibility test: for $L=10^5$ only $131\,065$ points pass it, $223$ of them give a perfect square, all $223$ pass the parity and distinctness tests, $155$ are smallest vertices, and all $155$ triangles have a perimeter below the limit. In our runs every version finished within ten seconds.

The total is a floating-point sum of $155$ positive terms below $10^5$, each made of three correctly rounded square roots of exactly represented integers. The worst-case error of such a sum grows only in proportion to the number of terms, ^[11] which keeps it below $10^{-6}$ here, while the exact total, recomputed with $50$-digit decimal arithmetic, lies more than $2\cdot10^{-5}$ away from the nearest rounding boundary, and a double-precision summation differs from it by less than $2\cdot10^{-9}$. The tolerance $10^{-12}$ in the perimeter test never matters, because no perimeter lies within $3700$ of $10^5$.

Besides the two checkpoints of the C++ program, the following was verified for this article:

A brute force that takes all triples of lattice points on each circle $x^2+y^2=n$ up to the radius allowed by Euler's formula and tests the orthocentre by its definition, through the two altitude conditions $(H-A)\cdot(B-C)=0$ and $(H-B)\cdot(A-C)=0$, reproduces the nine listed triangles and $291.0089$. For ten limits up to $1500$ it agrees triangle by triangle with the method, and the C++ program run with --perimeter= gives the same rounded sums; there are $13$ triangles with perimeter at most $100$ and $45$ with perimeter at most $1000$.
An independent full-size scan treats every lattice point of the disc of radius $r_{\max}+5$ as a vertex, finds the chord in floating point and confirms it with exact integers; at the vertex $H$ it takes every lattice diameter. It finds each of the $155$ triangles exactly once from each of its three vertices, and their largest circumradius is about $18521.35$.
For all $163$ triangles with perimeter at most $150\,000$, the side formula and Euler's formula of Section 8 hold exactly, and both perimeter bounds of Section 9 hold.
The C++ program built with Clang and with GCC, the Java program and the Python bridge print the same value, which is the correctly rounded exact total.

Footnotes and References

The references below supply the background facts used above. The reduction to a scan over one vertex, the conditions for a lattice chord, the radius bounds and all counts are derived explicitly in this article.

Project Euler 264 — Triangle Centres. The official statement considers the triangles with all vertices on lattice points, circumcentre at the origin O and orthocentre at H(5, 0), lists the nine of them with perimeter at most 50, whose perimeters add up to 291.0089 after rounding to four decimal places, and asks for the sum of the perimeters of all such triangles with perimeter at most 10⁵, rounded to four decimal places.
Wikipedia — Euler line. The article states Sylvester's relation OH = OA + OB + OC between the vectors from the circumcentre O to the orthocentre H and to the three vertices, and the formula OH² = 9R² − (a² + b² + c²) for the distance between the two centres.
Wikipedia — Thales's theorem. If A, B and C are distinct points on a circle and AC is a diameter, then the angle at B is a right angle.
Wikipedia — Parallelogram law. For all vectors x and y of a Euclidean space, 2‖x‖² + 2‖y‖² = ‖x + y‖² + ‖x − y‖².
Wikipedia — Euclid's lemma. In its general form, if an integer n divides a product ab and is coprime to a, then n divides b; the article proves this with Bézout's identity.
Wikipedia — Lexicographic order. On the Cartesian product of two ordered sets, pairs are compared by their first components and, when these are equal, by their second components; if both sets are totally ordered, the result is a total order.
Wikipedia — Concave function. Between any two points of its graph, the graph of a concave function of one variable lies on or above the straight line joining them; a twice differentiable function is concave exactly when its second derivative is never positive, and √x is an example.
C++ — Euler264.cpp. The linked immutable C++ revision contains lex_leq, isqrt_i64, dist, solve, round4 and run_checkpoints with the checkpoint 291.0089 for the limit 50. It is the implementation record for this article.
Wikipedia — Integer square root. The integer square root of a non-negative integer n is the greatest integer m with m² ≤ n, that is ⌊√n⌋.
Wikipedia — Long double. Most C compilers for x86 implement long double as the 80-bit extended precision type of the x86 hardware, while on 32-bit ARM and on 64-bit ARM under Windows and macOS long double has the same double-precision format as double.
Wikipedia — Kahan summation algorithm. Adding n floating-point numbers in sequence has a worst-case error that grows in proportion to n: the relative error bound of this naive summation is O(εn) times the condition number, which is 1 when all inputs are non-negative, and ε is about 10⁻¹⁶ in double precision.

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