Problem 998: Squaring the Triangle
View on Project EulerProject Euler Problem 998 Solution
EulerSolve provides an optimized solution for Project Euler Problem 998, Squaring the Triangle, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary For an integer-sided triangle \(\Delta\), let \(s(\Delta)\) be the side length of the smallest square that can contain \(\Delta\) after the triangle is rotated and translated in the plane. Project Euler 998 asks for the sum of the perimeters of all distinct integer-sided triangles whose minimum enclosing square has side length at most \(10^6\). The examples \(T(40)=346\), \(T(400)=76402\), and \(T(2000)=3237036\) are used as checkpoints. Mathematical Approach and Geometric Model Fix a square of side length \(L\), and rotate the whole configuration so the square is axis-aligned. A triangle that is tight for this square has vertices on the boundary. If one of its sides crosses the square with horizontal or vertical span \(L\), that side is the hypotenuse of a right triangle whose legs are \(L\) and some offset \(p\). Therefore its length is integral exactly when \[ h^2=L^2+p^2. \] This is the first major reduction: instead of searching arbitrary triangles, we enumerate Pythagorean segments attached to each possible square side \(L\). For every \(L\), the program stores pairs \((p,h)\), where \(0\le p\le L\) and \(h=\sqrt{L^2+p^2}\) is an integer. The restriction \(p\le L\) loses nothing because the complementary orientation is symmetric. Formal Derivation and Normal Form The enumeration rests on a normal-form statement....
Detailed mathematical approach
Problem Summary
For an integer-sided triangle \(\Delta\), let \(s(\Delta)\) be the side length of the smallest square that can contain \(\Delta\) after the triangle is rotated and translated in the plane. Project Euler 998 asks for the sum of the perimeters of all distinct integer-sided triangles whose minimum enclosing square has side length at most \(10^6\). The examples \(T(40)=346\), \(T(400)=76402\), and \(T(2000)=3237036\) are used as checkpoints.
Mathematical Approach and Geometric Model
Fix a square of side length \(L\), and rotate the whole configuration so the square is axis-aligned. A triangle that is tight for this square has vertices on the boundary. If one of its sides crosses the square with horizontal or vertical span \(L\), that side is the hypotenuse of a right triangle whose legs are \(L\) and some offset \(p\). Therefore its length is integral exactly when
\[ h^2=L^2+p^2. \]
This is the first major reduction: instead of searching arbitrary triangles, we enumerate Pythagorean segments attached to each possible square side \(L\). For every \(L\), the program stores pairs \((p,h)\), where \(0\le p\le L\) and \(h=\sqrt{L^2+p^2}\) is an integer. The restriction \(p\le L\) loses nothing because the complementary orientation is symmetric.
Formal Derivation and Normal Form
The enumeration rests on a normal-form statement. Let \(Q\) be a minimum square for a triangle. After rotating the plane, \(Q=[0,L]\times[0,L]\). A non-degenerate minimum cannot have all vertices strictly inside a smaller homothetic square, so at least two opposite supporting lines or two adjacent supporting lines are active. For a triangle, the supporting vertices can change only when a side direction becomes parallel to a square side. Thus every critical placement can be represented by boundary chords whose coordinate spans are determined by the side length \(L\) and by one boundary offset.
Lemma 1. If a side of a critical placement connects two boundary lines separated by \(L\), then its integral length is equivalent to a Pythagorean condition \(h^2=L^2+p^2\). This follows immediately from orthogonal projection: the side vector has components \(L\) and \(p\), up to sign. Conversely, every such Pythagorean segment can be placed on the square boundary and is therefore a legitimate building block for a candidate triangle.
Lemma 2. Pairing two boundary segments of the same square side \(L\) leaves only two combinatorial closures. Either the remaining side lies on a square boundary and has length \(p+q\), or it joins the two opposite residual endpoints and has squared length \((L-p)^2+(L-q)^2\). These are exactly the two families tested in the program. No third closure exists because the three vertices of the triangle occupy three chosen boundary endpoints, and after the first two side chords are fixed the last side is determined by which pair of residual endpoints is joined.
The inequality \(pq\ge L(L-p-q)\) is the algebraic form of the first closure's inside-square condition. It is obtained by writing the two boundary chords in coordinates, intersecting the two supporting rays from the boundary endpoints, and requiring the intersection coordinate to remain inside the square. This converts the geometric feasibility statement into one integer comparison, avoiding any floating-point decision in the first family.
Pythagorean Enumeration
The offsets are generated with Euclid's formula. Primitive triples are
\[ a=m^2-n^2,\qquad b=2mn,\qquad c=m^2+n^2, \]
with \(\gcd(m,n)=1\) and \(m-n\) odd. Each primitive triple is then scaled by \(k\). Either leg may play the role of the square side \(L\), and the other leg becomes the offset. Hence the code inserts both \((L, p, h)=(ka,kb,kc)\) and \((kb,ka,kc)\) whenever the offset is not larger than the side. Since the relevant leg is at most \(L\le 10^6\), the Euclid parameters only need to be tested up to \(\sqrt{2L}+2\). Duplicates are removed after sorting each list for a fixed \(L\).
First Family of Triangles
Choose two Pythagorean boundary segments for the same square side \(L\): \((p,h_p)\) and \((q,h_q)\). The first family is obtained when the remaining side of the triangle lies along one side of the square. Its length is
\[ b=p+q. \]
For the endpoints to lie on the boundary segment, we require \(b\le L\). The third vertex must also be on the correct side of the square rather than outside the \(L\times L\) box. In the boundary coordinates used by the implementation, this feasibility condition simplifies to
\[ pq\ge L(L-b). \]
When both inequalities hold, \((b,h_p,h_q)\) is a valid integer-sided triangle whose enclosing square has side length \(L\). Its perimeter is added after sorting the three side lengths into a canonical key.
Second Family of Triangles
The same two Pythagorean segments can also close through the opposite pair of endpoints. The remaining displacement has legs \(L-p\) and \(L-q\), so the third side is integral exactly when
\[ r^2=(L-p)^2+(L-q)^2. \]
This gives a candidate triangle \((h_p,h_q,r)\). Unlike the first family, the construction alone does not prove that \(L\) is the minimum enclosing square side. The same side triple might fit into a smaller square after a different rotation. Therefore every candidate from this family is passed to an independent minimum-square verifier, and it is accepted only when the verified value is at least \(L\), up to a small floating-point tolerance.
Minimum Enclosing Square Verification
For a triangle with side lengths \(a\le b\le c\), place the side \(c\) on the \(x\)-axis and compute the third vertex by the law of cosines:
\[ x=\frac{a^2+c^2-b^2}{2c},\qquad y=\sqrt{a^2-x^2}. \]
For a rotation angle \(\theta\), project the three vertices onto the rotated axes \(u\) and \(v\). The side length of the smallest axis-aligned square at that angle is
\[ w(\theta)=\max\{\max u-\min u,\ \max v-\min v\}. \]
The vertices that realize the maxima and minima change only when an edge becomes parallel to one of the square axes. Thus the interval \([0,\pi/2)\) is split at the edge directions. Inside each interval, both widths are fixed sinusoidal expressions. The optimum is found either at a boundary or at an interior angle where the two widths are equal. The routine minimum_square evaluates exactly these candidates, which is the standard rotating-calipers idea specialized to three points.
Deduplication and Summation
A triangle may be produced by several square sides, by both endpoint orderings, or by both geometric families. The problem asks for distinct triangles, so the program sorts every side triple and packs it into one integer key. The side lengths are below \(2^{21}\), so the key
\[ (a\ll 42)\,|\,(b\ll 21)\,|\,c \]
is collision-free for the target range. The perimeter is added only the first time a key appears.
Correctness Argument
Soundness. Every accepted triangle is integer-sided by construction. In the first family, the side lengths are \(p+q\), \(h_p\), and \(h_q\), with the Pythagorean equations ensuring integral hypotenuses and the feasibility inequality ensuring a placement inside the square. In the second family, the extra square test ensures that \(r\) is integral, and minimum_square verifies that the side triple does not belong to a smaller enclosing square. Therefore every inserted triangle satisfies \(s(\Delta)\le L\le n\).
Completeness. Take any integer-sided triangle with \(s(\Delta)\le n\), and place it in a minimum square of side \(L=s(\Delta)\). By the normal-form discussion, at least two of the tight side chords are Pythagorean boundary segments for the same \(L\). Their remaining endpoints close in one of the two possible combinatorial ways described in Lemma 2. Hence the triangle is generated by the corresponding iteration of the algorithm. If it appears through several placements, the canonical key keeps only one copy.
Numerical robustness. The enumeration, square tests, perimeter sums, and deduplication keys are exact integer operations. Floating point is used only in the secondary verifier for the second family. That verifier evaluates all critical angles of a three-point set and compares against \(L\) with a small tolerance; the published checkpoints provide an end-to-end guard against both missed candidates and false positives.
How the Code Works
pythagorean_offsets builds the list of all admissible Pythagorean boundary segments for every \(L\le n\). The main solve loop then considers all unordered pairs of segments attached to the same \(L\). It tests the first family with the algebraic inequality above, tests the second family with a square check and minimum_square, and inserts every accepted side triple into the global set. The three published checkpoints are asserted before the target value is printed.
Complexity Analysis
Let \(d_L\) be the number of Pythagorean offsets stored for side \(L\). The running time after triple generation is
\[ O\!\left(\sum_{L\le n} d_L^2\right), \]
because every unordered pair of offsets for the same \(L\) is tested once. The memory usage is \(O(\sum d_L+u)\), where \(u\) is the number of distinct accepted triangles. In practice the lists \(d_L\) are sparse, so \(n=10^6\) is feasible in compiled code.
References
- Problem page: Project Euler 998
- Wolfram MathWorld: Pythagorean Triple for Euclid's formula and primitive triples.
- Toussaint, Solving Geometric Problems with the Rotating Calipers for the support-line method used by the verifier.
- Wolfram MathWorld: Law of Cosines for reconstructing a triangle from its side lengths.
Problem 998 source code
C++
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <limits>
#include <numeric>
#include <tuple>
#include <unordered_set>
#include <utility>
#include <vector>
namespace {
using i64 = std::int64_t;
using u64 = std::uint64_t;
constexpr int TARGET = 1'000'000;
struct Leg {
int offset;
int hypotenuse;
};
i64 isqrt(const i64 n) {
i64 r = static_cast<i64>(std::sqrt(static_cast<long double>(n)));
while ((r + 1) * (r + 1) <= n) {
++r;
}
while (r * r > n) {
--r;
}
return r;
}
bool is_square(const i64 n, i64& root) {
root = isqrt(n);
return root * root == n;
}
long double normalize_angle(long double angle) {
const long double half_pi = std::acos(-1.0L) / 2.0L;
angle = std::fmod(angle, half_pi);
if (angle < 0) {
angle += half_pi;
}
if (angle < 1e-18L || half_pi - angle < 1e-18L) {
return 0.0L;
}
return angle;
}
struct WidthState {
long double width;
int u_min;
int u_max;
int v_min;
int v_max;
};
long double minimum_square(int a, int b, int c) {
if (a > b) {
std::swap(a, b);
}
if (b > c) {
std::swap(b, c);
}
if (a > b) {
std::swap(a, b);
}
const long double x = (static_cast<long double>(a) * a + static_cast<long double>(c) * c -
static_cast<long double>(b) * b) /
(2.0L * c);
long double y2 = static_cast<long double>(a) * a - x * x;
if (y2 < 0 && y2 > -1e-12L) {
y2 = 0;
}
const long double y = std::sqrt(y2);
const std::array<std::pair<long double, long double>, 3> points{{{0, 0}, {static_cast<long double>(c), 0}, {x, y}}};
auto evaluate = [&](const long double angle) {
const long double co = std::cos(angle);
const long double si = std::sin(angle);
std::array<long double, 3> u{};
std::array<long double, 3> v{};
for (int i = 0; i < 3; ++i) {
u[static_cast<std::size_t>(i)] = points[static_cast<std::size_t>(i)].first * co +
points[static_cast<std::size_t>(i)].second * si;
v[static_cast<std::size_t>(i)] = -points[static_cast<std::size_t>(i)].first * si +
points[static_cast<std::size_t>(i)].second * co;
}
const int u_min = static_cast<int>(std::min_element(u.begin(), u.end()) - u.begin());
const int u_max = static_cast<int>(std::max_element(u.begin(), u.end()) - u.begin());
const int v_min = static_cast<int>(std::min_element(v.begin(), v.end()) - v.begin());
const int v_max = static_cast<int>(std::max_element(v.begin(), v.end()) - v.begin());
return WidthState{std::max(u[static_cast<std::size_t>(u_max)] - u[static_cast<std::size_t>(u_min)],
v[static_cast<std::size_t>(v_max)] - v[static_cast<std::size_t>(v_min)]),
u_min,
u_max,
v_min,
v_max};
};
const long double half_pi = std::acos(-1.0L) / 2.0L;
std::vector<long double> angles{0};
for (int i = 0; i < 3; ++i) {
for (int j = i + 1; j < 3; ++j) {
const long double edge = std::atan2(points[static_cast<std::size_t>(j)].second -
points[static_cast<std::size_t>(i)].second,
points[static_cast<std::size_t>(j)].first -
points[static_cast<std::size_t>(i)].first);
angles.push_back(normalize_angle(edge));
angles.push_back(normalize_angle(edge + half_pi));
}
}
std::sort(angles.begin(), angles.end());
std::vector<long double> boundaries;
for (const long double angle : angles) {
if (boundaries.empty() || std::fabs(angle - boundaries.back()) > 1e-15L) {
boundaries.push_back(angle);
}
}
std::vector<long double> candidates = boundaries;
std::vector<long double> endpoints = boundaries;
endpoints.push_back(half_pi);
for (std::size_t i = 0; i + 1 < endpoints.size(); ++i) {
const long double lo = endpoints[i];
const long double hi = endpoints[i + 1];
if (hi - lo < 1e-15L) {
continue;
}
const WidthState state = evaluate((lo + hi) / 2.0L);
const auto& a0 = points[static_cast<std::size_t>(state.u_min)];
const auto& a1 = points[static_cast<std::size_t>(state.u_max)];
const auto& b0 = points[static_cast<std::size_t>(state.v_min)];
const auto& b1 = points[static_cast<std::size_t>(state.v_max)];
const long double ax = a1.first - a0.first;
const long double ay = a1.second - a0.second;
const long double bx = b1.first - b0.first;
const long double by = b1.second - b0.second;
const long double left = ax - by;
const long double right = ay + bx;
if (std::fabs(right) < 1e-18L) {
continue;
}
const long double root = std::atan2(-left, right);
for (const long double angle : {root, root + std::acos(-1.0L)}) {
const long double current = normalize_angle(angle);
if (lo + 1e-15L < current && current < hi - 1e-15L) {
candidates.push_back(current);
}
}
}
long double best = std::numeric_limits<long double>::max();
for (const long double angle : candidates) {
best = std::min(best, evaluate(angle).width);
}
return best;
}
u64 triangle_key(int a, int b, int c) {
if (a > b) {
std::swap(a, b);
}
if (b > c) {
std::swap(b, c);
}
if (a > b) {
std::swap(a, b);
}
return (static_cast<u64>(a) << 42) | (static_cast<u64>(b) << 21) | static_cast<u64>(c);
}
std::vector<std::vector<Leg>> pythagorean_offsets(const int limit) {
std::vector<std::vector<Leg>> by_side(static_cast<std::size_t>(limit + 1));
const int bound = isqrt(2LL * limit) + 2;
for (int m = 2; m <= bound; ++m) {
for (int n = 1; n < m; ++n) {
if (((m - n) & 1) == 0 || std::gcd(m, n) != 1) {
continue;
}
const int a = m * m - n * n;
const int b = 2 * m * n;
const int c = m * m + n * n;
for (const auto [side0, offset0] : {std::pair<int, int>{a, b}, {b, a}}) {
for (int k = 1; k * side0 <= limit; ++k) {
const int side = k * side0;
const int offset = k * offset0;
if (offset <= side) {
by_side[static_cast<std::size_t>(side)].push_back({offset, k * c});
}
}
}
}
}
for (std::vector<Leg>& legs : by_side) {
std::sort(legs.begin(), legs.end(), [](const Leg& lhs, const Leg& rhs) {
if (lhs.offset != rhs.offset) {
return lhs.offset < rhs.offset;
}
return lhs.hypotenuse < rhs.hypotenuse;
});
legs.erase(std::unique(legs.begin(), legs.end(), [](const Leg& lhs, const Leg& rhs) {
return lhs.offset == rhs.offset && lhs.hypotenuse == rhs.hypotenuse;
}),
legs.end());
}
return by_side;
}
u64 solve(const int limit) {
const std::vector<std::vector<Leg>> legs = pythagorean_offsets(limit);
std::unordered_set<u64> seen;
u64 total = 0;
auto add_triangle = [&](const int a, const int b, const int c) {
const u64 key = triangle_key(a, b, c);
if (seen.insert(key).second) {
total += static_cast<u64>(a + b + c);
}
};
for (int side = 1; side <= limit; ++side) {
const std::vector<Leg>& current = legs[static_cast<std::size_t>(side)];
for (std::size_t i = 0; i < current.size(); ++i) {
const i64 p = current[i].offset;
for (std::size_t j = i; j < current.size(); ++j) {
const i64 q = current[j].offset;
const i64 base = p + q;
if (base <= side && p * q >= static_cast<i64>(side) * (side - base)) {
add_triangle(static_cast<int>(base), current[i].hypotenuse, current[j].hypotenuse);
}
const i64 u = side - p;
const i64 v = side - q;
i64 third = 0;
if (u > 0 && v > 0 && is_square(u * u + v * v, third) &&
minimum_square(current[i].hypotenuse, current[j].hypotenuse, static_cast<int>(third)) + 1e-6L >=
side) {
add_triangle(current[i].hypotenuse, current[j].hypotenuse, static_cast<int>(third));
}
}
}
}
return total;
}
void run_checkpoints() {
assert(solve(40) == 346);
assert(solve(400) == 76'402);
assert(solve(2'000) == 3'237'036);
}
} // namespace
int main() {
run_checkpoints();
std::cout << solve(TARGET) << '\n';
return 0;
}
Python
from math import atan2, cos, fmod, gcd, isqrt, pi, sin, sqrt
TARGET = 1_000_000
HALF_PI = pi / 2.0
def is_square(n):
root = isqrt(n)
return root if root * root == n else -1
def normalize_angle(angle):
angle = fmod(angle, HALF_PI)
if angle < 0.0:
angle += HALF_PI
if angle < 1e-18 or HALF_PI - angle < 1e-18:
return 0.0
return angle
def minimum_square(a, b, c):
a, b, c = sorted((a, b, c))
x = (a * a + c * c - b * b) / (2.0 * c)
y2 = a * a - x * x
if -1e-12 < y2 < 0.0:
y2 = 0.0
y = sqrt(y2)
points = ((0.0, 0.0), (float(c), 0.0), (x, y))
def evaluate(angle):
co = cos(angle)
si = sin(angle)
u = [px * co + py * si for px, py in points]
v = [-px * si + py * co for px, py in points]
u_min = min(range(3), key=u.__getitem__)
u_max = max(range(3), key=u.__getitem__)
v_min = min(range(3), key=v.__getitem__)
v_max = max(range(3), key=v.__getitem__)
return (
max(u[u_max] - u[u_min], v[v_max] - v[v_min]),
u_min,
u_max,
v_min,
v_max,
)
angles = [0.0]
for i in range(3):
xi, yi = points[i]
for j in range(i + 1, 3):
xj, yj = points[j]
edge = atan2(yj - yi, xj - xi)
angles.append(normalize_angle(edge))
angles.append(normalize_angle(edge + HALF_PI))
angles.sort()
boundaries = []
for angle in angles:
if not boundaries or abs(angle - boundaries[-1]) > 1e-15:
boundaries.append(angle)
candidates = list(boundaries)
endpoints = boundaries + [HALF_PI]
for idx in range(len(endpoints) - 1):
lo = endpoints[idx]
hi = endpoints[idx + 1]
if hi - lo < 1e-15:
continue
_, u_min, u_max, v_min, v_max = evaluate((lo + hi) / 2.0)
ax = points[u_max][0] - points[u_min][0]
ay = points[u_max][1] - points[u_min][1]
bx = points[v_max][0] - points[v_min][0]
by = points[v_max][1] - points[v_min][1]
left = ax - by
right = ay + bx
if abs(right) < 1e-18:
continue
root = atan2(-left, right)
for angle in (root, root + pi):
current = normalize_angle(angle)
if lo + 1e-15 < current < hi - 1e-15:
candidates.append(current)
return min(evaluate(angle)[0] for angle in candidates)
def triangle_key(a, b, c):
a, b, c = sorted((a, b, c))
return (a << 42) | (b << 21) | c
def pythagorean_offsets(limit):
by_side = [[] for _ in range(limit + 1)]
bound = isqrt(2 * limit) + 2
for m in range(2, bound + 1):
for n in range(1, m):
if ((m - n) & 1) == 0 or gcd(m, n) != 1:
continue
a = m * m - n * n
b = 2 * m * n
c = m * m + n * n
for side0, offset0 in ((a, b), (b, a)):
side = side0
offset = offset0
hypotenuse = c
while side <= limit:
if offset <= side:
by_side[side].append((offset, hypotenuse))
side += side0
offset += offset0
hypotenuse += c
for side in range(limit + 1):
if by_side[side]:
by_side[side] = sorted(set(by_side[side]))
return by_side
def solve(limit):
legs = pythagorean_offsets(limit)
seen = set()
total = 0
def add_triangle(a, b, c):
nonlocal total
key = triangle_key(a, b, c)
if key not in seen:
seen.add(key)
total += a + b + c
for side, current in enumerate(legs):
size = len(current)
for i in range(size):
p, hyp_i = current[i]
for j in range(i, size):
q, hyp_j = current[j]
base = p + q
if base <= side and p * q >= side * (side - base):
add_triangle(base, hyp_i, hyp_j)
u = side - p
v = side - q
if u > 0 and v > 0:
third = is_square(u * u + v * v)
if third >= 0 and minimum_square(hyp_i, hyp_j, third) + 1e-6 >= side:
add_triangle(hyp_i, hyp_j, third)
return total
def run_checkpoints():
assert solve(40) == 346
assert solve(400) == 76_402
assert solve(2_000) == 3_237_036
if __name__ == "__main__":
run_checkpoints()
print(solve(TARGET))
Java
import java.util.ArrayList;
import java.util.Collections;
import java.util.Comparator;
import java.util.HashSet;
import java.util.List;
import java.util.Set;
public class Euler998 {
private static final int TARGET = 1_000_000;
private static final double HALF_PI = Math.PI / 2.0;
private static final class Leg {
final int offset;
final int hypotenuse;
Leg(int offset, int hypotenuse) {
this.offset = offset;
this.hypotenuse = hypotenuse;
}
}
private static final class WidthState {
final double width;
final int uMin;
final int uMax;
final int vMin;
final int vMax;
WidthState(double width, int uMin, int uMax, int vMin, int vMax) {
this.width = width;
this.uMin = uMin;
this.uMax = uMax;
this.vMin = vMin;
this.vMax = vMax;
}
}
private static long isqrt(long n) {
long r = (long) Math.sqrt((double) n);
while ((r + 1) * (r + 1) <= n) {
++r;
}
while (r * r > n) {
--r;
}
return r;
}
private static long isSquare(long n) {
long root = isqrt(n);
return root * root == n ? root : -1;
}
private static int gcd(int a, int b) {
while (b != 0) {
int t = a % b;
a = b;
b = t;
}
return Math.abs(a);
}
private static double normalizeAngle(double angle) {
angle = angle % HALF_PI;
if (angle < 0.0) {
angle += HALF_PI;
}
if (angle < 1e-18 || HALF_PI - angle < 1e-18) {
return 0.0;
}
return angle;
}
private static WidthState evaluate(double[][] points, double angle) {
double co = Math.cos(angle);
double si = Math.sin(angle);
double[] u = new double[3];
double[] v = new double[3];
for (int i = 0; i < 3; ++i) {
u[i] = points[i][0] * co + points[i][1] * si;
v[i] = -points[i][0] * si + points[i][1] * co;
}
int uMin = 0;
int uMax = 0;
int vMin = 0;
int vMax = 0;
for (int i = 1; i < 3; ++i) {
if (u[i] < u[uMin]) {
uMin = i;
}
if (u[i] > u[uMax]) {
uMax = i;
}
if (v[i] < v[vMin]) {
vMin = i;
}
if (v[i] > v[vMax]) {
vMax = i;
}
}
double width = Math.max(u[uMax] - u[uMin], v[vMax] - v[vMin]);
return new WidthState(width, uMin, uMax, vMin, vMax);
}
private static double minimumSquare(int a, int b, int c) {
if (a > b) {
int t = a;
a = b;
b = t;
}
if (b > c) {
int t = b;
b = c;
c = t;
}
if (a > b) {
int t = a;
a = b;
b = t;
}
double x = ((double) a * a + (double) c * c - (double) b * b) / (2.0 * c);
double y2 = (double) a * a - x * x;
if (y2 < 0.0 && y2 > -1e-12) {
y2 = 0.0;
}
double y = Math.sqrt(y2);
double[][] points = {
{0.0, 0.0},
{(double) c, 0.0},
{x, y}
};
List<Double> angles = new ArrayList<>();
angles.add(0.0);
for (int i = 0; i < 3; ++i) {
for (int j = i + 1; j < 3; ++j) {
double edge = Math.atan2(points[j][1] - points[i][1], points[j][0] - points[i][0]);
angles.add(normalizeAngle(edge));
angles.add(normalizeAngle(edge + HALF_PI));
}
}
Collections.sort(angles);
List<Double> boundaries = new ArrayList<>();
for (double angle : angles) {
if (boundaries.isEmpty() || Math.abs(angle - boundaries.get(boundaries.size() - 1)) > 1e-15) {
boundaries.add(angle);
}
}
List<Double> candidates = new ArrayList<>(boundaries);
List<Double> endpoints = new ArrayList<>(boundaries);
endpoints.add(HALF_PI);
for (int i = 0; i + 1 < endpoints.size(); ++i) {
double lo = endpoints.get(i);
double hi = endpoints.get(i + 1);
if (hi - lo < 1e-15) {
continue;
}
WidthState state = evaluate(points, (lo + hi) / 2.0);
double ax = points[state.uMax][0] - points[state.uMin][0];
double ay = points[state.uMax][1] - points[state.uMin][1];
double bx = points[state.vMax][0] - points[state.vMin][0];
double by = points[state.vMax][1] - points[state.vMin][1];
double left = ax - by;
double right = ay + bx;
if (Math.abs(right) < 1e-18) {
continue;
}
double root = Math.atan2(-left, right);
for (double angle : new double[] {root, root + Math.PI}) {
double current = normalizeAngle(angle);
if (lo + 1e-15 < current && current < hi - 1e-15) {
candidates.add(current);
}
}
}
double best = Double.MAX_VALUE;
for (double angle : candidates) {
best = Math.min(best, evaluate(points, angle).width);
}
return best;
}
private static long triangleKey(int a, int b, int c) {
if (a > b) {
int t = a;
a = b;
b = t;
}
if (b > c) {
int t = b;
b = c;
c = t;
}
if (a > b) {
int t = a;
a = b;
b = t;
}
return ((long) a << 42) | ((long) b << 21) | (long) c;
}
@SuppressWarnings("unchecked")
private static List<Leg>[] pythagoreanOffsets(int limit) {
List<Leg>[] bySide = new ArrayList[limit + 1];
for (int i = 0; i <= limit; ++i) {
bySide[i] = new ArrayList<>();
}
int bound = (int) isqrt(2L * limit) + 2;
for (int m = 2; m <= bound; ++m) {
for (int n = 1; n < m; ++n) {
if (((m - n) & 1) == 0 || gcd(m, n) != 1) {
continue;
}
int a = m * m - n * n;
int b = 2 * m * n;
int c = m * m + n * n;
int[] sides = {a, b};
int[] offsets = {b, a};
for (int idx = 0; idx < 2; ++idx) {
int side0 = sides[idx];
int offset0 = offsets[idx];
int side = side0;
int offset = offset0;
int hypotenuse = c;
while (side <= limit) {
if (offset <= side) {
bySide[side].add(new Leg(offset, hypotenuse));
}
side += side0;
offset += offset0;
hypotenuse += c;
}
}
}
}
Comparator<Leg> cmp = Comparator.comparingInt((Leg leg) -> leg.offset)
.thenComparingInt(leg -> leg.hypotenuse);
for (int side = 0; side <= limit; ++side) {
List<Leg> legs = bySide[side];
if (legs.isEmpty()) {
continue;
}
legs.sort(cmp);
List<Leg> unique = new ArrayList<>();
Leg prev = null;
for (Leg leg : legs) {
if (prev == null || prev.offset != leg.offset || prev.hypotenuse != leg.hypotenuse) {
unique.add(leg);
prev = leg;
}
}
bySide[side] = unique;
}
return bySide;
}
private static long solve(int limit) {
List<Leg>[] legs = pythagoreanOffsets(limit);
Set<Long> seen = new HashSet<>();
long total = 0;
for (int side = 1; side <= limit; ++side) {
List<Leg> current = legs[side];
for (int i = 0; i < current.size(); ++i) {
long p = current.get(i).offset;
int hypI = current.get(i).hypotenuse;
for (int j = i; j < current.size(); ++j) {
long q = current.get(j).offset;
int hypJ = current.get(j).hypotenuse;
long base = p + q;
if (base <= side && p * q >= (long) side * (side - base)) {
long key = triangleKey((int) base, hypI, hypJ);
if (seen.add(key)) {
total += base + hypI + hypJ;
}
}
long u = side - p;
long v = side - q;
if (u > 0 && v > 0) {
long third = isSquare(u * u + v * v);
if (third >= 0 && minimumSquare(hypI, hypJ, (int) third) + 1e-6 >= side) {
long key = triangleKey(hypI, hypJ, (int) third);
if (seen.add(key)) {
total += hypI + hypJ + third;
}
}
}
}
}
}
return total;
}
private static void check(long actual, long expected, String label) {
if (actual != expected) {
throw new IllegalStateException(label + ": expected " + expected + ", got " + actual);
}
}
private static void runCheckpoints() {
check(solve(40), 346L, "T(40)");
check(solve(400), 76_402L, "T(400)");
check(solve(2_000), 3_237_036L, "T(2000)");
}
public static void main(String[] args) {
runCheckpoints();
System.out.println(solve(TARGET));
}
}