Problem 562: Maximal Perimeter

Project Euler Problem 562 Solution

EulerSolve provides an optimized solution for Project Euler Problem 562, Maximal Perimeter, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We study integer-coordinate triangles $A,B,C$ whose vertices lie in the closed disk $$x^2+y^2\le r^2.$$ Among all such triangles with area $$\Delta=\frac12,$$ the task is to maximize the perimeter. If several triangles reach the same perimeter, the one with larger circumradius $R$ is preferred, and the reported quantity is $$T(r)=\frac{R}{r}.$$ The brute-force search space is enormous, so the solution rewrites the geometry in terms of a primitive base edge and a one-dimensional lattice family for the third vertex. Mathematical Approach Fix one candidate triangle and choose the side $AB$ as the base. Write $$u=B-A=(a,b),\qquad w=C-A,\qquad L=|u|=\sqrt{a^2+b^2}.$$ The implementations only keep configurations where the chosen base is the longest side, so $$|w|\le L,\qquad |u-w|\le L.$$ Step 1: Convert the area condition into a determinant equation For lattice points, the doubled area is the absolute determinant: $$2\Delta=|\det(u,w)|.$$ Because $\Delta=\frac12$, every admissible triangle satisfies $$|\det(u,w)|=1.$$ If $A=(x_A,y_A)$ and $C=(x,y)$, then $w=C-A$, so the same condition becomes $$a y-b x=a y_A-b x_A\pm1.$$ This immediately forces $\gcd(a,b)=1$; if $a$ and $b$ had a common factor $g>1$, then the left-hand side would always be divisible by $g$, so it could never equal $\pm1$....

Detailed mathematical approach

Problem Summary

We study integer-coordinate triangles $A,B,C$ whose vertices lie in the closed disk

$$x^2+y^2\le r^2.$$

Among all such triangles with area

$$\Delta=\frac12,$$

the task is to maximize the perimeter. If several triangles reach the same perimeter, the one with larger circumradius $R$ is preferred, and the reported quantity is

$$T(r)=\frac{R}{r}.$$

The brute-force search space is enormous, so the solution rewrites the geometry in terms of a primitive base edge and a one-dimensional lattice family for the third vertex.

Mathematical Approach

Fix one candidate triangle and choose the side $AB$ as the base. Write

$$u=B-A=(a,b),\qquad w=C-A,\qquad L=|u|=\sqrt{a^2+b^2}.$$

The implementations only keep configurations where the chosen base is the longest side, so

$$|w|\le L,\qquad |u-w|\le L.$$

Step 1: Convert the area condition into a determinant equation

For lattice points, the doubled area is the absolute determinant:

$$2\Delta=|\det(u,w)|.$$

Because $\Delta=\frac12$, every admissible triangle satisfies

$$|\det(u,w)|=1.$$

If $A=(x_A,y_A)$ and $C=(x,y)$, then $w=C-A$, so the same condition becomes

$$a y-b x=a y_A-b x_A\pm1.$$

This immediately forces $\gcd(a,b)=1$; if $a$ and $b$ had a common factor $g>1$, then the left-hand side would always be divisible by $g$, so it could never equal $\pm1$.

Step 2: Parameterize all admissible third vertices

Since $\gcd(a,b)=1$, the linear Diophantine equation

$$-b x+a y=1$$

has an integer solution. If the required right-hand side is some integer $k$, one particular point on the line is obtained by scaling that solution by $k$. Every other integer point on the same line differs by a multiple of the base vector $u$, so the full family is

$$C(t)=C_0+t(a,b),\qquad t\in\mathbb{Z}.$$

Thus, once $A$, $B$, and the sign choice $\pm1$ are fixed, the third vertex can only move along one lattice line parallel to the base.

Step 3: Intersect the lattice line with the disk

Substituting $C(t)=C_0+t u$ into the disk constraint gives

$$|C_0+t u|^2\le r^2,$$

which expands to the quadratic inequality

$$L^2 t^2+2(C_0\cdot u)t+\bigl(|C_0|^2-r^2\bigr)\le0.$$

The discriminant tells us whether that line meets the disk at all. When it does, the admissible integers form one contiguous interval

$$t_{\min}\le t\le t_{\max}.$$

The implementations compute those endpoints with exact integer floor and ceiling arithmetic, so the feasibility test does not depend on floating-point rounding.

Step 4: Express the two other sides through one integer dot product

Let

$$q=u\cdot w.$$

This is an integer because both vectors have integer coordinates. Lagrange's identity gives

$$|u|^2|w|^2=(u\cdot w)^2+\det(u,w)^2=q^2+1,$$

hence

$$|w|=\frac{\sqrt{q^2+1}}{L}.$$

For the third side, use $u-w=B-C$. Then

$$u\cdot(u-w)=L^2-q,\qquad \det(u,u-w)=\det(u,w),$$

$$|u-w|=\frac{\sqrt{(L^2-q)^2+1}}{L}.$$

The perimeter therefore becomes

$$P=L+\frac{\sqrt{q^2+1}+\sqrt{(L^2-q)^2+1}}{L}.$$

Step 5: Derive the pruning bound $P\le2L+\frac1L$

Because the chosen base must be the longest side, both $q$ and $L^2-q$ are positive integers. For every integer $n\ge1$,

$$\sqrt{n^2+1}\le n+\frac{1}{2n}.$$

Apply this inequality once to $q$ and once to $L^2-q$:

$$\sqrt{q^2+1}+\sqrt{(L^2-q)^2+1} \le q+\frac{1}{2q}+(L^2-q)+\frac{1}{2(L^2-q)} \le L^2+1.$$

After dividing by $L$, we obtain

$$|w|+|u-w|\le L+\frac1L,$$

and therefore

$$\boxed{P\le2L+\frac1L.}$$

This is the key pruning inequality used by the search. Once even this optimistic upper bound drops below the best perimeter already found, no shorter base can improve the answer.

Step 6: Search edges by deficiency from the diameter

Any edge inside the disk has length at most $2r$. For that reason, promising base edges are those with $L$ very close to $2r$. The implementations measure this gap by the deficiency

$$d=4r^2-L^2.$$

Small $d$ means an almost diametral base, so candidate primitive edges are generated in increasing deficiency, or equivalently decreasing $L$.

If the current best perimeter is $P^*$, then a future edge can only win if

$$2L+\frac1L>P^*.$$

Solving the equality gives the larger root

$$L_{\min}=\frac{P^*+\sqrt{(P^*)^2-8}}{4}.$$

Hence no edge with $L<L_{\min}$ can matter, and the search only needs to cover deficiencies up to

$$d_{\max}=4r^2-L_{\min}^2.$$

Worked Example

Take

$$A=(0,0),\qquad B=(2,1).$$

Then $u=(2,1)$ and $L^2=5$. The area condition becomes

$$|\det(u,C-A)|=|2y-x|=1.$$

Choose the $+1$ branch. One integer point on the line $2y-x=1$ is $(-1,0)$, so every solution has the form

$$C(t)=(-1,0)+t(2,1).$$

Inside the disk of radius $3$, we need

$$(-1+2t)^2+t^2\le9,$$

which leaves $t=0$ and $t=1$. The choice $t=1$ gives $C=(1,1)$, and then

$$|C-A|^2=2,\qquad |C-B|^2=1,\qquad P=\sqrt5+\sqrt2+1.$$

Since $\Delta=\frac12$, the circumradius is

$$R=\frac{|AB|\cdot|AC|\cdot|BC|}{4\Delta} =\frac{\sqrt5\cdot\sqrt2\cdot1}{2} =\frac{\sqrt{10}}{2}.$$

How the Code Works

The implementation first precomputes, for each integer $x\in[-r,r]$, the largest $|y|$ still inside the disk. That allows fast testing of whether both endpoints of a translated base edge stay inside the circle.

Next it builds primitive edge vectors $(a,b)$ with $1\le b\le a$, so rotational and reflection symmetries are not revisited. Only vectors whose squared length lies in the current deficiency window are kept, and they are processed from longer to shorter edges.

For each admissible translation of the base, the C++, Python, and Java implementations examine the two determinant lines corresponding to area $+\frac12$ and $-\frac12$. An integer point on the chosen line is produced with the extended Euclidean algorithm, then shifted along the base direction so the quadratic disk test stays numerically small and exact.

Every integer point in the resulting $t$-interval is checked. The implementation discards cases where the third point coincides with an endpoint or where one of the other two sides exceeds the base length. Surviving candidates are compared first by perimeter and then, on ties, by circumradius.

The C++ version evaluates edge batches in parallel, the Java version mirrors the same exact search with arbitrary-precision arithmetic where needed, and the Python version acts as a thin bridge that launches the native computation and extracts the final numeric answer.

Complexity Analysis

Let $E(D)$ be the number of primitive edge vectors whose deficiency is at most $D$. Precomputing the disk boundary table costs $O(r)$ time and $O(r)$ memory. Generating and sorting the current candidate edges costs $O(E(D)\log E(D))$.

For a single edge, the solver scans all translations that keep both endpoints in the disk; in the worst case this can be as large as $O(r^2)$. Each valid translation yields at most two determinant lines, and each line contributes one contiguous interval of integer parameters for the third vertex. A deliberately loose bound for one search round is therefore $O(r^2E(D))$, but in practice the deficiency window stays small, the perimeter bound removes many edges early, and the integer $t$-interval is usually short.

Overall, the method is strongly data-dependent but memory usage remains modest: $O(r+E(D))$, plus thread-local state in the parallel C++ version.

Footnotes and References

Problem page: https://projecteuler.net/problem=562
Pick's theorem: Wikipedia - Pick's theorem
Circumradius formulas: Wikipedia - Circumscribed circle
Extended Euclidean algorithm: Wikipedia - Extended Euclidean algorithm
Linear Diophantine equations: Wikipedia - Diophantine equation

Problem 562 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <limits>
#include <numeric>
#include <string>
#include <thread>
#include <vector>

namespace {

using i64 = std::int64_t;
using u64 = std::uint64_t;
using i128 = __int128_t;
using u128 = __uint128_t;

constexpr long double kPerimeterEps = 1e-15L;

struct TriangleResult {
    long double perimeter = -1.0L;
    long double circumradius = -1.0L;
    i64 ax = 0;
    i64 ay = 0;
    i64 bx = 0;
    i64 by = 0;
    i64 cx = 0;
    i64 cy = 0;
};

struct Edge {
    i64 a = 0;
    i64 b = 0;
    u64 s2 = 0;
};

i128 floor_div(i128 a, i128 b) {
    i128 q = a / b;
    i128 r = a % b;
    if (r != 0 && ((r > 0) != (b > 0))) --q;
    return q;
}

i128 ceil_div(i128 a, i128 b) {
    i128 q = a / b;
    i128 r = a % b;
    if (r != 0 && ((r > 0) == (b > 0))) ++q;
    return q;
}

i128 round_div_nearest(i128 a, i128 b) {
    // Rounds a/b to nearest integer, ties away from zero.
    if (a >= 0) {
        return (a + b / 2) / b;
    }
    return -((-a + b / 2) / b);
}

u64 isqrt_u64(u64 n) {
    long double x = std::sqrt(static_cast<long double>(n));
    u64 r = static_cast<u64>(x);
    while ((r + 1) * (r + 1) <= n) ++r;
    while (r * r > n) --r;
    return r;
}

u64 ceil_sqrt_u64(u64 n) {
    u64 r = isqrt_u64(n);
    if (r * r < n) ++r;
    return r;
}

u128 isqrt_u128(u128 n) {
    if (n == 0) return 0;
    long double approx = std::sqrt(static_cast<long double>(n));
    u128 x = static_cast<u128>(approx);
    if (x == 0) x = 1;
    while ((x + 1) * (x + 1) <= n) ++x;
    while (x * x > n) --x;
    return x;
}

i64 extended_gcd(i64 a, i64 b, i64 &x, i64 &y) {
    // Returns g=gcd(a,b), and x,y such that a*x + b*y = g.
    i64 x0 = 1, y0 = 0;
    i64 x1 = 0, y1 = 1;
    i64 aa = a, bb = b;
    while (bb != 0) {
        i64 q = aa / bb;
        i64 t = aa - q * bb;
        aa = bb;
        bb = t;

        i64 nx = x0 - q * x1;
        i64 ny = y0 - q * y1;
        x0 = x1;
        y0 = y1;
        x1 = nx;
        y1 = ny;
    }
    x = x0;
    y = y0;
    return aa;
}

long double perimeter_upper_bound(u64 s2) {
    long double L = std::sqrt(static_cast<long double>(s2));
    return 2.0L * L + 1.0L / L;
}

bool better_triangle(const TriangleResult &cand, const TriangleResult &best) {
    if (cand.perimeter > best.perimeter + kPerimeterEps) return true;
    if (std::fabsl(cand.perimeter - best.perimeter) <= kPerimeterEps &&
        cand.circumradius > best.circumradius + 1e-18L) {
        return true;
    }
    return false;
}

class Euler562Solver {
public:
    Euler562Solver(i64 r, int threads)
        : r_(r), r2_(static_cast<u64>(r) * static_cast<u64>(r)),
          threads_(std::max(1, threads)) {
        precompute_ymax();
    }

    TriangleResult solve_exact() {
        TriangleResult best;
        u64 d_bound = 128;
        const u64 max_d = 4ULL * r2_;

        while (true) {
            if (d_bound > max_d) d_bound = max_d;
            TriangleResult cur = search_with_deficiency_bound(d_bound, best);
            if (better_triangle(cur, best)) best = cur;

            if (best.perimeter < 0.0L) {
                if (d_bound == max_d) break;
                d_bound = std::min<u64>(max_d, d_bound * 2ULL);
                continue;
            }

            const long double p = best.perimeter;
            const long double disc = p * p - 8.0L;
            const long double l_bound = (p + std::sqrt(disc)) / 4.0L;
            const long double d_req_ld =
                4.0L * static_cast<long double>(r2_) - l_bound * l_bound;
            u64 d_req = 0;
            if (d_req_ld > 0.0L) {
                d_req = static_cast<u64>(std::ceill(d_req_ld - 1e-18L));
            }

            if (d_bound >= d_req || d_bound == max_d) break;
            d_bound = std::max<u64>(std::min<u64>(max_d, d_bound * 2ULL), d_req);
        }

        return best;
    }

private:
    i64 r_;
    u64 r2_;
    int threads_;
    std::vector<int> y_max_;

    void precompute_ymax() {
        y_max_.assign(static_cast<std::size_t>(2 * r_ + 1), 0);
        for (i64 x = -r_; x <= r_; ++x) {
            u64 rem = r2_ - static_cast<u64>(x * x);
            y_max_[static_cast<std::size_t>(x + r_)] = static_cast<int>(isqrt_u64(rem));
        }
    }

    TriangleResult search_with_deficiency_bound(u64 d_bound,
                                                const TriangleResult &initial_best) const {
        std::vector<Edge> edges;
        build_candidate_edges(d_bound, initial_best.perimeter, edges);
        if (edges.empty()) return initial_best;

        std::sort(edges.begin(), edges.end(), [](const Edge &lhs, const Edge &rhs) {
            if (lhs.s2 != rhs.s2) return lhs.s2 > rhs.s2;
            if (lhs.a != rhs.a) return lhs.a < rhs.a;
            return lhs.b < rhs.b;
        });

        TriangleResult best = initial_best;

        const std::size_t batch_size = 2048;
        std::size_t i = 0;
        while (i < edges.size()) {
            if (best.perimeter >= 0.0L &&
                perimeter_upper_bound(edges[i].s2) <= best.perimeter + kPerimeterEps) {
                break;
            }

            std::size_t j = std::min(edges.size(), i + batch_size);
            const int use_threads =
                (threads_ <= 1 || (j - i) < 64) ? 1 : std::min<int>(threads_, static_cast<int>(j - i));

            if (use_threads == 1) {
                for (std::size_t idx = i; idx < j; ++idx) {
                    const Edge &e = edges[idx];
                    if (best.perimeter >= 0.0L &&
                        perimeter_upper_bound(e.s2) <= best.perimeter + kPerimeterEps) {
                        continue;
                    }
                    TriangleResult cand = evaluate_edge(e, best.perimeter);
                    if (better_triangle(cand, best)) best = cand;
                }
            } else {
                std::vector<TriangleResult> locals(static_cast<std::size_t>(use_threads), best);
                std::vector<std::thread> pool;
                pool.reserve(static_cast<std::size_t>(use_threads));

                for (int t = 0; t < use_threads; ++t) {
                    pool.emplace_back([&, t]() {
                        TriangleResult local_best = best;
                        for (std::size_t idx = i + static_cast<std::size_t>(t); idx < j;
                             idx += static_cast<std::size_t>(use_threads)) {
                            const Edge &e = edges[idx];
                            if (local_best.perimeter >= 0.0L &&
                                perimeter_upper_bound(e.s2) <= local_best.perimeter + kPerimeterEps) {
                                continue;
                            }
                            TriangleResult cand = evaluate_edge(e, local_best.perimeter);
                            if (better_triangle(cand, local_best)) local_best = cand;
                        }
                        locals[static_cast<std::size_t>(t)] = local_best;
                    });
                }
                for (auto &th : pool) th.join();
                for (const auto &cand : locals) {
                    if (better_triangle(cand, best)) best = cand;
                }
            }

            i = j;
        }

        return best;
    }

    void build_candidate_edges(u64 d_bound, long double best_perimeter,
                               std::vector<Edge> &edges) const {
        edges.clear();
        const u64 max_s2 = 4ULL * r2_;
        const u64 min_s2 = (d_bound >= max_s2) ? 0ULL : (max_s2 - d_bound);

        // By symmetry we only need 1 <= b <= a.
        const u64 a_min_sq = (min_s2 + 1ULL) / 2ULL;
        i64 a_min = static_cast<i64>(ceil_sqrt_u64(a_min_sq));
        if (a_min < 1) a_min = 1;

        for (i64 a = a_min; a <= 2 * r_; ++a) {
            const u64 aa = static_cast<u64>(a) * static_cast<u64>(a);
            if (aa > max_s2) break;

            const u64 rem_hi = max_s2 - aa;
            i64 b_max = static_cast<i64>(isqrt_u64(rem_hi));
            if (b_max > a) b_max = a;
            if (b_max < 1) continue;

            i64 b_min = 1;
            if (aa < min_s2) {
                const u64 rem_lo = min_s2 - aa;
                b_min = static_cast<i64>(ceil_sqrt_u64(rem_lo));
            }
            if (b_min > b_max) continue;

            for (i64 b = b_min; b <= b_max; ++b) {
                if (std::gcd(a, b) != 1) continue;
                u64 s2 = aa + static_cast<u64>(b) * static_cast<u64>(b);
                if (s2 < min_s2 || s2 > max_s2) continue;
                if (best_perimeter >= 0.0L &&
                    perimeter_upper_bound(s2) <= best_perimeter + kPerimeterEps) {
                    continue;
                }
                edges.push_back({a, b, s2});
            }
        }
    }

    TriangleResult evaluate_edge(const Edge &e, long double threshold) const {
        TriangleResult best;
        best.perimeter = threshold;
        best.circumradius = -1.0L;

        const i64 a = e.a;
        const i64 b = e.b;
        const u64 s2_u64 = e.s2;
        const i128 s2 = static_cast<i128>(s2_u64);
        const long double L = std::sqrt(static_cast<long double>(s2_u64));

        if (perimeter_upper_bound(s2_u64) <= threshold + kPerimeterEps) return best;

        i64 x1 = 0, y1 = 0;
        i64 g = extended_gcd(-b, a, x1, y1);
        if (g == -1) {
            x1 = -x1;
            y1 = -y1;
        } else if (g != 1) {
            return best;
        }

        const i64 x_lo = std::max(-r_, -r_ - a);
        const i64 x_hi = std::min(r_, r_ - a);
        if (x_lo > x_hi) return best;

        for (i64 xA = x_lo; xA <= x_hi; ++xA) {
            const i64 xB = xA + a;
            const int y1_max = y_max_[static_cast<std::size_t>(xA + r_)];
            const int y2_max = y_max_[static_cast<std::size_t>(xB + r_)];
            const i64 y_lo = std::max<i64>(-y1_max, -static_cast<i64>(y2_max) - b);
            const i64 y_hi = std::min<i64>(y1_max, static_cast<i64>(y2_max) - b);
            if (y_lo > y_hi) continue;

            for (i64 yA = y_lo; yA <= y_hi; ++yA) {
                const i64 yB = yA + b;
                const i128 detA = static_cast<i128>(a) * yA - static_cast<i128>(b) * xA;

                for (int sgn : {-1, 1}) {
                    const i128 k = detA + static_cast<i128>(sgn);

                    // Raw point on line det(u, C)=k: C_raw = k * (x1, y1).
                    i128 x_raw = static_cast<i128>(x1) * k;
                    i128 y_raw = static_cast<i128>(y1) * k;

                    // Shift along u to nearest point to origin to keep numbers bounded.
                    i128 dot = x_raw * a + y_raw * b;
                    i128 shift = round_div_nearest(-dot, s2);
                    i128 x0 = x_raw + shift * a;
                    i128 y0 = y_raw + shift * b;

                    const i128 Acoef = s2;
                    const i128 Bcoef = 2 * (x0 * a + y0 * b);
                    const i128 Ccoef = x0 * x0 + y0 * y0 - static_cast<i128>(r2_);

                    i128 disc = Bcoef * Bcoef - 4 * Acoef * Ccoef;
                    if (disc < 0) continue;

                    u128 sqrt_disc = isqrt_u128(static_cast<u128>(disc));
                    i128 denom = 2 * Acoef;
                    i128 t_min = ceil_div(-Bcoef - static_cast<i128>(sqrt_disc), denom);
                    i128 t_max = floor_div(-Bcoef + static_cast<i128>(sqrt_disc), denom);
                    if (t_min > t_max) continue;

                    for (i128 t = t_min; t <= t_max; ++t) {
                        i128 xC128 = x0 + t * a;
                        i128 yC128 = y0 + t * b;
                        i128 normC2 = xC128 * xC128 + yC128 * yC128;
                        if (normC2 > static_cast<i128>(r2_)) continue;

                        i64 xC = static_cast<i64>(xC128);
                        i64 yC = static_cast<i64>(yC128);

                        if ((xC == xA && yC == yA) || (xC == xB && yC == yB)) continue;

                        i128 dx1 = xC128 - static_cast<i128>(xA);
                        i128 dy1 = yC128 - static_cast<i128>(yA);
                        i128 d2_128 = dx1 * dx1 + dy1 * dy1;

                        i128 dx2 = xC128 - static_cast<i128>(xB);
                        i128 dy2 = yC128 - static_cast<i128>(yB);
                        i128 d3_128 = dx2 * dx2 + dy2 * dy2;

                        if (d2_128 > s2 || d3_128 > s2) continue;  // u is longest side.

                        u64 d2 = static_cast<u64>(d2_128);
                        u64 d3 = static_cast<u64>(d3_128);

                        long double p =
                            L + std::sqrt(static_cast<long double>(d2)) + std::sqrt(static_cast<long double>(d3));
                        if (p <= best.perimeter + kPerimeterEps) continue;

                        long double R =
                            std::sqrt(static_cast<long double>(s2_u64) * static_cast<long double>(d2) *
                                      static_cast<long double>(d3)) /
                            2.0L;

                        TriangleResult cand;
                        cand.perimeter = p;
                        cand.circumradius = R;
                        cand.ax = xA;
                        cand.ay = yA;
                        cand.bx = xB;
                        cand.by = yB;
                        cand.cx = xC;
                        cand.cy = yC;

                        if (better_triangle(cand, best)) best = cand;
                    }
                }
            }
        }

        return best;
    }
};

long double compute_T(i64 r, int threads, TriangleResult *out_triangle = nullptr) {
    Euler562Solver solver(r, threads);
    TriangleResult best = solver.solve_exact();
    if (out_triangle) *out_triangle = best;
    return best.circumradius / static_cast<long double>(r);
}

bool run_validation(int threads) {
    struct Checkpoint {
        i64 r;
        long double expected;
        long double tol;
    };

    const std::vector<Checkpoint> checks = {
        {5, std::sqrt(19669.0L / 50.0L), 1e-12L},
        {10, 97.26729L, 1e-5L},
        {100, 9157.64707L, 1e-5L},
    };

    for (const auto &cp : checks) {
        TriangleResult tri;
        long double got = compute_T(cp.r, threads, &tri);
        long double err = std::fabsl(got - cp.expected);
        std::cout << "T(" << cp.r << ") = " << std::setprecision(12) << std::fixed << got << "\n";
        if (err > cp.tol) {
            std::cerr << "Validation failed at r=" << cp.r << ": expected ~" << cp.expected
                      << ", got " << got << "\n";
            std::cerr << "Triangle: A=(" << tri.ax << "," << tri.ay << "), B=(" << tri.bx << ","
                      << tri.by << "), C=(" << tri.cx << "," << tri.cy << ")\n";
            return false;
        }
    }

    std::cerr << "Validation checkpoints passed.\n";
    return true;
}

}  // namespace

int main(int argc, char **argv) {
    i64 target_r = 10'000'000;
    int threads = static_cast<int>(std::thread::hardware_concurrency());
    if (threads <= 0) threads = 1;
    bool validate = true;

    for (int i = 1; i < argc; ++i) {
        std::string arg = argv[i];
        if ((arg == "--r" || arg == "-r") && i + 1 < argc) {
            target_r = std::stoll(argv[++i]);
        } else if ((arg == "--threads" || arg == "-t") && i + 1 < argc) {
            threads = std::max(1, std::stoi(argv[++i]));
        } else if (arg == "--no-validation") {
            validate = false;
        }
    }

    if (validate) {
        if (!run_validation(std::min(threads, 4))) return 1;
    }

    TriangleResult tri;
    long double T = compute_T(target_r, threads, &tri);
    unsigned long long rounded = static_cast<unsigned long long>(std::llround(T));

    std::cout << std::setprecision(12) << std::fixed;
    std::cout << "T(" << target_r << ") = " << T << "\n";
    std::cout << "Rounded answer = " << rounded << "\n";
    std::cout << rounded << "\n";

    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

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


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


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


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

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


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

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

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


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

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.Collections;
import java.util.List;

public class Euler562 {
    static final double kPerimeterEps = 1e-15;

    static class TriangleResult {
        double perimeter = -1.0;
        double circumradius = -1.0;
        long ax, ay, bx, by, cx, cy;
    }

    static class Edge implements Comparable<Edge> {
        long a, b;
        long s2;

        Edge(long a, long b, long s2) {
            this.a = a;
            this.b = b;
            this.s2 = s2;
        }

        @Override
        public int compareTo(Edge o) {
            if (this.s2 != o.s2)
                return this.s2 > o.s2 ? -1 : 1;
            if (this.a != o.a)
                return this.a < o.a ? -1 : 1;
            return this.b < o.b ? -1 : (this.b == o.b ? 0 : 1);
        }
    }

    static BigInteger floorDiv(BigInteger a, BigInteger b) {
        BigInteger[] qr = a.divideAndRemainder(b);
        BigInteger q = qr[0];
        BigInteger r = qr[1];
        if (r.signum() != 0 && ((r.signum() > 0) != (b.signum() > 0))) {
            q = q.subtract(BigInteger.ONE);
        }
        return q;
    }

    static BigInteger ceilDiv(BigInteger a, BigInteger b) {
        BigInteger[] qr = a.divideAndRemainder(b);
        BigInteger q = qr[0];
        BigInteger r = qr[1];
        if (r.signum() != 0 && ((r.signum() > 0) == (b.signum() > 0))) {
            q = q.add(BigInteger.ONE);
        }
        return q;
    }

    static BigInteger roundDivNearest(BigInteger a, BigInteger b) {
        if (a.signum() >= 0) {
            return a.add(b.divide(BigInteger.valueOf(2))).divide(b);
        }
        return a.negate().add(b.divide(BigInteger.valueOf(2))).divide(b).negate();
    }

    static long isqrt(long n) {
        if (n < 0)
            return 0;
        long r = (long) Math.sqrt(n);
        while ((r + 1) * (r + 1) <= n)
            r++;
        while (r * r > n)
            r--;
        return r;
    }

    static long ceilSqrt(long n) {
        long r = isqrt(n);
        if (r * r < n)
            r++;
        return r;
    }

    static BigInteger isqrtBig(BigInteger n) {
        if (n.signum() == 0)
            return BigInteger.ZERO;
        BigInteger r = BigInteger.valueOf((long) Math.sqrt(n.doubleValue()));
        if (r.signum() == 0)
            r = BigInteger.ONE;
        while (r.add(BigInteger.ONE).multiply(r.add(BigInteger.ONE)).compareTo(n) <= 0)
            r = r.add(BigInteger.ONE);
        while (r.multiply(r).compareTo(n) > 0)
            r = r.subtract(BigInteger.ONE);
        return r;
    }

    static long gcd(long a, long b) {
        a = Math.abs(a);
        b = Math.abs(b);
        while (b != 0) {
            long temp = b;
            b = a % b;
            a = temp;
        }
        return a;
    }

    static long[] extendedGcd(long a, long b) {
        long x0 = 1, y0 = 0;
        long x1 = 0, y1 = 1;
        long aa = a, bb = b;
        while (bb != 0) {
            long q = aa / bb;
            long t = aa - q * bb;
            aa = bb;
            bb = t;

            long nx = x0 - q * x1;
            long ny = y0 - q * y1;
            x0 = x1;
            y0 = y1;
            x1 = nx;
            y1 = ny;
        }
        return new long[] { x0, y0, aa };
    }

    static double perimeterUpperBound(long s2) {
        double L = Math.sqrt((double) s2);
        return 2.0 * L + 1.0 / L;
    }

    static boolean betterTriangle(TriangleResult cand, TriangleResult best) {
        if (cand.perimeter > best.perimeter + kPerimeterEps)
            return true;
        if (Math.abs(cand.perimeter - best.perimeter) <= kPerimeterEps && cand.circumradius > best.circumradius + 1e-18)
            return true;
        return false;
    }

    static class Solver {
        long r;
        long r2;
        int[] yMax;

        Solver(long r) {
            this.r = r;
            this.r2 = r * r;
            yMax = new int[(int) (2 * r + 1)];
            for (long x = -r; x <= r; x++) {
                long rem = r2 - x * x;
                yMax[(int) (x + r)] = (int) isqrt(rem);
            }
        }

        TriangleResult solveExact() {
            TriangleResult best = new TriangleResult();
            long dBound = 128;
            long maxD = 4 * r2;

            while (true) {
                if (dBound > maxD)
                    dBound = maxD;
                TriangleResult cur = searchWithDeficiencyBound(dBound, best);
                if (betterTriangle(cur, best))
                    best = cur;

                if (best.perimeter < 0.0) {
                    if (dBound == maxD)
                        break;
                    dBound = Math.min(maxD, dBound * 2);
                    continue;
                }

                double p = best.perimeter;
                double disc = p * p - 8.0;
                double lBound = (p + Math.sqrt(disc)) / 4.0;
                double dReqLd = 4.0 * r2 - lBound * lBound;
                long dReq = 0;
                if (dReqLd > 0.0) {
                    dReq = (long) Math.ceil(dReqLd - 1e-18);
                }

                if (dBound >= dReq || dBound == maxD)
                    break;
                dBound = Math.max(Math.min(maxD, dBound * 2), dReq);
            }

            return best;
        }

        TriangleResult searchWithDeficiencyBound(long dBound, TriangleResult initialBest) {
            List<Edge> edges = new ArrayList<>();
            buildCandidateEdges(dBound, initialBest.perimeter, edges);
            if (edges.isEmpty())
                return initialBest;

            Collections.sort(edges);

            TriangleResult best = initialBest;

            for (Edge e : edges) {
                if (best.perimeter >= 0.0 && perimeterUpperBound(e.s2) <= best.perimeter + kPerimeterEps) {
                    break;
                }
                TriangleResult cand = evaluateEdge(e, best.perimeter);
                if (betterTriangle(cand, best))
                    best = cand;
            }

            return best;
        }

        void buildCandidateEdges(long dBound, double bestPerimeter, List<Edge> edges) {
            long maxS2 = 4 * r2;
            long minS2 = (dBound >= maxS2) ? 0 : (maxS2 - dBound);

            long aMinSq = (minS2 + 1) / 2;
            long aMin = ceilSqrt(aMinSq);
            if (aMin < 1)
                aMin = 1;

            for (long a = aMin; a <= 2 * r; a++) {
                long aa = a * a;
                if (aa > maxS2)
                    break;

                long remHi = maxS2 - aa;
                long bMax = Math.min(a, isqrt(remHi));
                if (bMax < 1)
                    continue;

                long bMin = 1;
                if (aa < minS2) {
                    bMin = ceilSqrt(minS2 - aa);
                }
                if (bMin > bMax)
                    continue;

                for (long b = bMin; b <= bMax; b++) {
                    if (gcd(a, b) != 1)
                        continue;
                    long s2 = aa + b * b;
                    if (s2 < minS2 || s2 > maxS2)
                        continue;
                    if (bestPerimeter >= 0.0 && perimeterUpperBound(s2) <= bestPerimeter + kPerimeterEps) {
                        continue;
                    }
                    edges.add(new Edge(a, b, s2));
                }
            }
        }

        TriangleResult evaluateEdge(Edge e, double threshold) {
            TriangleResult best = new TriangleResult();
            best.perimeter = threshold;

            long a = e.a;
            long b = e.b;
            long s2 = e.s2;
            BigInteger s2Big = BigInteger.valueOf(s2);
            double L = Math.sqrt(s2);

            if (perimeterUpperBound(s2) <= threshold + kPerimeterEps)
                return best;

            long[] extGcd = extendedGcd(-b, a);
            long x1 = extGcd[0];
            long y1 = extGcd[1];
            long g = extGcd[2];

            if (g == -1) {
                x1 = -x1;
                y1 = -y1;
            } else if (g != 1) {
                return best;
            }

            long xLo = Math.max(-r, -r - a);
            long xHi = Math.min(r, r - a);
            if (xLo > xHi)
                return best;

            for (long xA = xLo; xA <= xHi; xA++) {
                long xB = xA + a;
                int y1Max = yMax[(int) (xA + r)];
                int y2Max = yMax[(int) (xB + r)];
                long yLo = Math.max(-y1Max, -y2Max - b);
                long yHi = Math.min(y1Max, y2Max - b);
                if (yLo > yHi)
                    continue;

                for (long yA = yLo; yA <= yHi; yA++) {
                    long yB = yA + b;
                    BigInteger detA = BigInteger.valueOf(a).multiply(BigInteger.valueOf(yA))
                            .subtract(BigInteger.valueOf(b).multiply(BigInteger.valueOf(xA)));

                    for (int sgn : new int[] { -1, 1 }) {
                        BigInteger k = detA.add(BigInteger.valueOf(sgn));

                        BigInteger xRaw = BigInteger.valueOf(x1).multiply(k);
                        BigInteger yRaw = BigInteger.valueOf(y1).multiply(k);

                        BigInteger dot = xRaw.multiply(BigInteger.valueOf(a)).add(yRaw.multiply(BigInteger.valueOf(b)));
                        BigInteger shift = roundDivNearest(dot.negate(), s2Big);
                        BigInteger x0 = xRaw.add(shift.multiply(BigInteger.valueOf(a)));
                        BigInteger y0 = yRaw.add(shift.multiply(BigInteger.valueOf(b)));

                        BigInteger Acoef = s2Big;
                        BigInteger Bcoef = x0.multiply(BigInteger.valueOf(a)).add(y0.multiply(BigInteger.valueOf(b)))
                                .multiply(BigInteger.valueOf(2));
                        BigInteger Ccoef = x0.multiply(x0).add(y0.multiply(y0)).subtract(BigInteger.valueOf(r2));

                        BigInteger disc = Bcoef.multiply(Bcoef)
                                .subtract(Acoef.multiply(Ccoef).multiply(BigInteger.valueOf(4)));
                        if (disc.signum() < 0)
                            continue;

                        BigInteger sqrtDisc = isqrtBig(disc);
                        BigInteger denom = Acoef.multiply(BigInteger.valueOf(2));

                        BigInteger tMin = floorDiv(Bcoef.negate().subtract(sqrtDisc), denom).add(BigInteger.ONE);
                        BigInteger tMax = floorDiv(Bcoef.negate().add(sqrtDisc), denom);
                        if (Bcoef.negate().subtract(sqrtDisc).remainder(denom).signum() == 0)
                            tMin = tMin.subtract(BigInteger.ONE);

                        tMin = ceilDiv(Bcoef.negate().subtract(sqrtDisc), denom);

                        if (tMin.compareTo(tMax) > 0)
                            continue;

                        // Because the loop logic works perfectly, iterate inside:
                        long minT = tMin.longValue();
                        long maxT = tMax.longValue();

                        for (long t = minT; t <= maxT; t++) {
                            BigInteger tBig = BigInteger.valueOf(t);
                            BigInteger xC = x0.add(tBig.multiply(BigInteger.valueOf(a)));
                            BigInteger yC = y0.add(tBig.multiply(BigInteger.valueOf(b)));

                            if (xC.multiply(xC).add(yC.multiply(yC)).compareTo(BigInteger.valueOf(r2)) > 0)
                                continue;

                            long xC_L = xC.longValue();
                            long yC_L = yC.longValue();

                            if ((xC_L == xA && yC_L == yA) || (xC_L == xB && yC_L == yB))
                                continue;

                            long dx1 = xC_L - xA;
                            long dy1 = yC_L - yA;
                            long d2 = dx1 * dx1 + dy1 * dy1;

                            long dx2 = xC_L - xB;
                            long dy2 = yC_L - yB;
                            long d3 = dx2 * dx2 + dy2 * dy2;

                            if (d2 > s2 || d3 > s2)
                                continue;

                            double p = L + Math.sqrt(d2) + Math.sqrt(d3);
                            if (p <= best.perimeter + kPerimeterEps)
                                continue;

                            double R = Math.sqrt((double) s2 * d2 * d3) / 2.0;

                            TriangleResult cand = new TriangleResult();
                            cand.perimeter = p;
                            cand.circumradius = R;
                            if (betterTriangle(cand, best)) {
                                best = cand;
                            }
                        }
                    }
                }
            }
            return best;
        }
    }

    public static String solve() {
        Solver solver = new Solver(10000000);
        TriangleResult best = solver.solveExact();
        double T = best.circumradius / 10000000.0;
        return Long.toString(Math.round(T));
    }

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