Problem 385: Ellipses Inside Triangles
View on Project EulerProject Euler Problem 385 Solution
EulerSolve provides an optimized solution for Project Euler Problem 385, Ellipses Inside Triangles, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary The local C++ validator shows that every admissible triangle can be represented in centered form $$P_1=(x_1,y_1),\qquad P_2=(x_2,y_2),\qquad P_3=(-x_1-x_2,-y_1-y_2),$$ so automatically $$x_1+x_2+x_3=0,\qquad y_1+y_2+y_3=0.$$ The ellipse condition is already converted by the implementation into the Diophantine system $$ (2x_1+x_2)y_1 + (x_1+2x_2)y_2 = 0, $$ $$ x_1^2+x_1x_2+x_2^2 - (y_1^2+y_1y_2+y_2^2) = 39. $$ We must sum the Euclidean areas of all non-degenerate integer triangles whose coordinates lie inside \([-N,N]^2\), for the real target \(N=10^9\). The brute-force search in solve_bruteforce_small is only feasible for tiny \(N\), so the fast solver rewrites the problem as a finite family of Pell-type recurrences. Mathematical Approach 1. Extract a primitive direction from the \(x\)-coordinates Write the three \(x\)-coordinates as $$x_1=du,\qquad x_2=dv,\qquad x_3=-d(u+v),$$ where \(d>0\) and \(\gcd(u,v)=1\). The linear constraint then becomes $$a y_1 + b y_2 = 0,\qquad a=2u+v,\qquad b=u+2v.$$ All integer solutions of this linear equation have the form $$y_1=\frac{b}{m}k,\qquad y_2=-\frac{a}{m}k,\qquad y_3=-y_1-y_2=\frac{u-v}{m}k,$$ where \(m=\gcd(a,b)\) and \(k\in\mathbb{Z}\). Because $$2a-b=3u,\qquad 2b-a=3v,$$ every common divisor of \(a\) and \(b\) must divide \(3\). Since \((u,v)\) is primitive, \(m\in\{1,3\}\)....
Detailed mathematical approach
Problem Summary
The local C++ validator shows that every admissible triangle can be represented in centered form
$$P_1=(x_1,y_1),\qquad P_2=(x_2,y_2),\qquad P_3=(-x_1-x_2,-y_1-y_2),$$
so automatically
$$x_1+x_2+x_3=0,\qquad y_1+y_2+y_3=0.$$
The ellipse condition is already converted by the implementation into the Diophantine system
$$ (2x_1+x_2)y_1 + (x_1+2x_2)y_2 = 0, $$
$$ x_1^2+x_1x_2+x_2^2 - (y_1^2+y_1y_2+y_2^2) = 39. $$
We must sum the Euclidean areas of all non-degenerate integer triangles whose coordinates lie inside \([-N,N]^2\), for the real target \(N=10^9\). The brute-force search in solve_bruteforce_small is only feasible for tiny \(N\), so the fast solver rewrites the problem as a finite family of Pell-type recurrences.
Mathematical Approach
1. Extract a primitive direction from the \(x\)-coordinates
Write the three \(x\)-coordinates as
$$x_1=du,\qquad x_2=dv,\qquad x_3=-d(u+v),$$
where \(d>0\) and \(\gcd(u,v)=1\). The linear constraint then becomes
$$a y_1 + b y_2 = 0,\qquad a=2u+v,\qquad b=u+2v.$$
All integer solutions of this linear equation have the form
$$y_1=\frac{b}{m}k,\qquad y_2=-\frac{a}{m}k,\qquad y_3=-y_1-y_2=\frac{u-v}{m}k,$$
where \(m=\gcd(a,b)\) and \(k\in\mathbb{Z}\). Because
$$2a-b=3u,\qquad 2b-a=3v,$$
every common divisor of \(a\) and \(b\) must divide \(3\). Since \((u,v)\) is primitive, \(m\in\{1,3\}\). The admissible primitive directions kept by the program always satisfy \(m=3\).
2. Substitute into the quadratic constraint
Define the norm
$$q_0=u^2+uv+v^2.$$
Then the \(x\)-part of the quadratic equation becomes
$$x_1^2+x_1x_2+x_2^2=d^2q_0.$$
For the \(y\)-part, using \(m=3\), \(a=2u+v\), and \(b=u+2v\), we obtain
$$y_1^2+y_1y_2+y_2^2=\frac{b^2-ab+a^2}{9}k^2.$$
Now
$$a^2-ab+b^2=3(u^2+uv+v^2)=3q_0,$$
so
$$y_1^2+y_1y_2+y_2^2=\frac{q_0}{3}k^2.$$
Substituting into the validator equation gives
$$39=q_0d^2-\frac{q_0}{3}k^2=\frac{q_0}{3}(3d^2-k^2),$$
hence
$$q_0(3d^2-k^2)=117,\qquad k^2-3d^2=-\frac{117}{q_0}.$$
This identity is the algebraic backbone of the fast method.
3. Why only \(q_0=3\) and \(q_0=39\)
The previous formula shows that \(q_0\) must divide \(117\). The code then enumerates primitive pairs \((u,v)\) satisfying
$$u^2+uv+v^2=q_0,\qquad \gcd(u,v)=1,\qquad \gcd(2u+v,u+2v)=3.$$
A direct enumeration leaves exactly two norm families:
$$q_0=3 \quad\text{with 6 primitive directions},\qquad q_0=39 \quad\text{with 12 primitive directions}.$$
No admissible primitive direction exists for \(q_0=117\), so the whole search collapses to \(18\) direction templates. That is why build_directions loops only over q0 in (3, 39).
4. Pell equations and boundary cutoffs
The two surviving norms lead to
$$q_0=3 \Longrightarrow k^2-3d^2=-39,\qquad q_0=39 \Longrightarrow k^2-3d^2=-3.$$
The implementations use the smallest positive seeds found by inspection:
$$ (d,k)=(4,3),(5,6)\ \text{for }-39,\qquad (d,k)=(2,3)\ \text{for }-3. $$
All larger solutions are generated by multiplication with the fundamental unit \(2+\sqrt{3}\), which yields the recurrence
$$k_{n+1}=2k_n+3d_n,\qquad d_{n+1}=k_n+2d_n.$$
For one fixed direction \((u,v)\), the coordinate bounds are
$$x_{\mathrm{scale}}=\max(|u|,|v|,|u+v|),$$
$$y_{\mathrm{scale}}=\frac{\max(|2u+v|,|u+2v|,|u-v|)}{3},$$
so the code keeps only pairs with
$$d\le \left\lfloor\frac{N}{x_{\mathrm{scale}}}\right\rfloor,\qquad k\le \left\lfloor\frac{N}{y_{\mathrm{scale}}}\right\rfloor.$$
5. Area formula, worked example, and deduplication
Using the explicit coordinates, the doubled area simplifies to
$$2A=\left|x_1(y_2-y_3)+x_2(y_3-y_1)+x_3(y_1-y_2)\right|=2q_0|dk|,$$
so the actual area is
$$\boxed{A=q_0|dk|}.$$
Take the primitive direction \((u,v)=(1,1)\), which belongs to the \(q_0=3\) family. With the seed \((d,k)=(4,3)\), we get
$$ (x_1,y_1)=(4,3),\qquad (x_2,y_2)=(4,-3),\qquad (x_3,y_3)=(-8,0), $$
and therefore
$$A=3\cdot 4\cdot 3=36.$$
This is one of the triangles contributing to the checkpoint value \(72\) at \(N=8\). Because different direction/sign choices can reconstruct the same triangle, the implementations sort the three vertices and use the sorted triple as a hash key. If a duplicate is seen, the code checks that the area is identical before merging.
How the Code Works
The three solution files implement the same pipeline. build_directions enumerates the \(18\) primitive direction templates and computes their personal \((d_{\max},k_{\max})\) cutoffs. generate_pell_pairs precomputes all Pell solutions for the two negative right-hand sides up to the largest cutoff required by any direction. build_triangle_key reconstructs the three vertices, rejects out-of-range or degenerate cases, and returns the sorted vertex triple together with the closed-form area \(q_0|dk|\).
The main loop selects the correct Pell list for each direction, tries both signs of \(k\), and inserts the triangle into a hash map keyed by its sorted vertices. The C++ version adds extra validation: checkpoint values \(72\), \(252\), \(34632\), and \(3529008\) for \(N=8,10,100,1000\), plus a brute-force cross-check against solve_bruteforce_small(20). The Python and Java versions use the same fast construction without the checkpoint harness.
Complexity Analysis
For this problem the number of direction templates is constant: \(18\). Each Pell sequence grows exponentially because the recurrence multiplies by \(2+\sqrt{3}\), so the number of generated \((d,k)\) pairs below the bounds is \(O(\log N)\) per family. The total running time is therefore proportional to the number of generated candidates and hash insertions, rather than to any polynomial scan over lattice coordinates. Memory usage is \(O(T)\), where \(T\) is the number of distinct triangles stored in the hash map.
Footnotes and References
- Problem page: https://projecteuler.net/problem=385
- Primary derivation source in this repository:
solutionsCpp/Euler385.cpp,solutionsPython/Euler385.py,solutionsJava/Euler385.java - Pell equation: Wikipedia — Pell equation
- Binary quadratic forms: Wikipedia — Binary quadratic form
- Triangle area / determinant formula: Wikipedia — Shoelace formula
Problem 385 source code
C++
#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <limits>
#include <numeric>
#include <string>
#include <thread>
#include <unordered_map>
#include <utility>
#include <vector>
namespace {
using i64 = std::int64_t;
using u64 = std::uint64_t;
using u128 = unsigned __int128;
constexpr i64 kDefaultN = 1'000'000'000LL;
struct Options {
i64 n = kDefaultN;
bool run_checkpoints = true;
bool allow_multithreading = true;
unsigned requested_threads = 0U;
};
struct Point {
i64 x = 0;
i64 y = 0;
bool operator<(const Point& other) const {
if (x != other.x) {
return x < other.x;
}
return y < other.y;
}
bool operator==(const Point& other) const {
return x == other.x && y == other.y;
}
};
struct TriangleKey {
std::array<Point, 3> vertices{};
bool operator==(const TriangleKey& other) const {
return vertices == other.vertices;
}
};
struct TriangleKeyHash {
std::size_t operator()(const TriangleKey& key) const {
std::size_t h = 0ULL;
for (const Point& p : key.vertices) {
const std::size_t hx = std::hash<i64>{}(p.x);
const std::size_t hy = std::hash<i64>{}(p.y);
h ^= hx + 0x9e3779b97f4a7c15ULL + (h << 6U) + (h >> 2U);
h ^= hy + 0x9e3779b97f4a7c15ULL + (h << 6U) + (h >> 2U);
}
return h;
}
};
struct Direction {
int q0 = 0;
int u = 0;
int v = 0;
int m = 3;
int pell_n = 0; // k^2 - 3 d^2 = pell_n.
i64 d_max = 0;
i64 k_max = 0;
};
bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& value) {
const std::string p(prefix);
if (arg.rfind(p, 0) != 0U) {
return false;
}
const std::string tail = arg.substr(p.size());
if (tail.empty()) {
return false;
}
u64 parsed = 0ULL;
for (const char c : tail) {
if (c < '0' || c > '9') {
return false;
}
const u64 digit = static_cast<u64>(c - '0');
if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
return false;
}
parsed = parsed * 10ULL + digit;
}
value = parsed;
return true;
}
bool parse_unsigned_after_prefix(const std::string& arg,
const char* prefix,
unsigned& value) {
u64 parsed = 0ULL;
if (!parse_u64_after_prefix(arg, prefix, parsed)) {
return false;
}
if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
return false;
}
value = static_cast<unsigned>(parsed);
return true;
}
bool parse_arguments(const int argc, char** argv, Options& options) {
for (int i = 1; i < argc; ++i) {
const std::string arg(argv[i]);
if (arg == "--skip-checkpoints") {
options.run_checkpoints = false;
continue;
}
if (arg == "--single-thread") {
options.allow_multithreading = false;
continue;
}
u64 parsed_u64 = 0ULL;
if (parse_u64_after_prefix(arg, "--n=", parsed_u64)) {
if (parsed_u64 > static_cast<u64>(std::numeric_limits<i64>::max())) {
std::cerr << "--n is too large for this implementation.\n";
return false;
}
options.n = static_cast<i64>(parsed_u64);
continue;
}
unsigned parsed_unsigned = 0U;
if (parse_unsigned_after_prefix(arg, "--threads=", parsed_unsigned)) {
options.requested_threads = parsed_unsigned;
continue;
}
std::cerr << "Unknown argument: " << arg << '\n';
return false;
}
return true;
}
unsigned choose_thread_count(const bool allow_multithreading,
const unsigned requested_threads,
const std::size_t workload_units) {
constexpr std::size_t kMinUnitsForParallel = 8ULL;
if (!allow_multithreading || workload_units < 2ULL || workload_units < kMinUnitsForParallel) {
return 1U;
}
unsigned threads = requested_threads;
if (threads == 0U) {
threads = std::thread::hardware_concurrency();
if (threads == 0U) {
threads = 1U;
}
}
return std::max(1U, std::min<unsigned>(threads, static_cast<unsigned>(workload_units)));
}
i64 abs_i64(const i64 x) {
return x >= 0 ? x : -x;
}
u64 abs_to_u64(const i64 x) {
return static_cast<u64>(x >= 0 ? x : -x);
}
std::string to_string_u128(u128 value) {
if (value == 0U) {
return "0";
}
std::string digits;
while (value > 0U) {
const unsigned digit = static_cast<unsigned>(value % 10U);
digits.push_back(static_cast<char>('0' + digit));
value /= 10U;
}
std::reverse(digits.begin(), digits.end());
return digits;
}
std::vector<std::pair<i64, i64>> generate_pell_pairs(const int pell_n,
const i64 d_max,
const i64 k_max) {
// Problem-specific seeds for k^2 - 3 d^2 = -39 and = -3.
std::vector<std::pair<i64, i64>> seeds;
if (pell_n == -39) {
seeds = {{3, 4}, {6, 5}};
} else if (pell_n == -3) {
seeds = {{3, 2}};
} else {
return {};
}
std::vector<std::pair<i64, i64>> all;
for (const auto [k0, d0] : seeds) {
i64 k = k0;
i64 d = d0;
while (d <= d_max && k <= k_max) {
all.emplace_back(d, k);
if (k > (std::numeric_limits<i64>::max() - 3LL * d) / 2LL ||
d > (std::numeric_limits<i64>::max() - k) / 2LL) {
break;
}
const i64 next_k = 2LL * k + 3LL * d;
const i64 next_d = k + 2LL * d;
k = next_k;
d = next_d;
}
}
std::sort(all.begin(), all.end());
all.erase(std::unique(all.begin(), all.end()), all.end());
return all;
}
std::vector<Direction> build_directions(const i64 n) {
std::vector<Direction> directions;
for (const int q0 : {3, 39}) {
const int pell_n = -117 / q0; // Derived from q0(3 d^2 - k^2) = 117.
const int limit = static_cast<int>(std::sqrt(static_cast<long double>(q0))) + 2;
for (int u = -limit; u <= limit; ++u) {
for (int v = -limit; v <= limit; ++v) {
if (u == 0 && v == 0) {
continue;
}
if (std::gcd(abs_i64(u), abs_i64(v)) != 1) {
continue;
}
if (u * u + u * v + v * v != q0) {
continue;
}
const int a = 2 * u + v;
const int b = u + 2 * v;
const int m = static_cast<int>(std::gcd(abs_i64(a), abs_i64(b)));
if (m != 3) {
continue;
}
const i64 x_scale = std::max({abs_i64(u), abs_i64(v), abs_i64(u + v)});
const i64 y_scale_raw = std::max({abs_i64(a), abs_i64(b), abs_i64(u - v)});
if (x_scale == 0 || y_scale_raw % m != 0) {
continue;
}
const i64 y_scale = y_scale_raw / m;
if (y_scale == 0) {
continue;
}
Direction dir;
dir.q0 = q0;
dir.u = u;
dir.v = v;
dir.m = m;
dir.pell_n = pell_n;
dir.d_max = n / x_scale;
dir.k_max = n / y_scale;
if (dir.d_max > 0 && dir.k_max > 0) {
directions.push_back(dir);
}
}
}
}
std::sort(directions.begin(),
directions.end(),
[](const Direction& lhs, const Direction& rhs) {
if (lhs.q0 != rhs.q0) {
return lhs.q0 < rhs.q0;
}
if (lhs.u != rhs.u) {
return lhs.u < rhs.u;
}
return lhs.v < rhs.v;
});
return directions;
}
bool build_triangle_key(const Direction& dir,
const i64 d,
const i64 k,
const i64 n,
TriangleKey& out_key,
u64& out_area) {
const i64 u = dir.u;
const i64 v = dir.v;
const i64 m = dir.m;
const i64 a = 2LL * u + v;
const i64 b = u + 2LL * v;
if (a % m != 0LL || b % m != 0LL || (u - v) % m != 0LL) {
return false;
}
const i64 x1 = d * u;
const i64 x2 = d * v;
const i64 x3 = -d * (u + v);
const i64 y1 = (b / m) * k;
const i64 y2 = -(a / m) * k;
const i64 y3 = ((u - v) / m) * k;
if (abs_i64(x1) > n || abs_i64(x2) > n || abs_i64(x3) > n || abs_i64(y1) > n ||
abs_i64(y2) > n || abs_i64(y3) > n) {
return false;
}
std::array<Point, 3> points = {{{x1, y1}, {x2, y2}, {x3, y3}}};
std::sort(points.begin(), points.end());
if (points[0] == points[1] || points[1] == points[2]) {
return false;
}
out_key.vertices = points;
out_area = static_cast<u64>(dir.q0) * abs_to_u64(d * k);
return true;
}
u128 sum_triangle_areas_from_map(const std::unordered_map<TriangleKey, u64, TriangleKeyHash>& map) {
u128 total = 0U;
for (const auto& entry : map) {
total += static_cast<u128>(entry.second);
}
return total;
}
u128 solve_fast(const i64 n,
const bool allow_multithreading,
const unsigned requested_threads,
std::size_t* out_triangle_count = nullptr) {
const std::vector<Direction> directions = build_directions(n);
if (directions.empty()) {
if (out_triangle_count != nullptr) {
*out_triangle_count = 0ULL;
}
return 0U;
}
i64 max_d_n39 = 0;
i64 max_k_n39 = 0;
i64 max_d_n3 = 0;
i64 max_k_n3 = 0;
for (const Direction& dir : directions) {
if (dir.pell_n == -39) {
max_d_n39 = std::max(max_d_n39, dir.d_max);
max_k_n39 = std::max(max_k_n39, dir.k_max);
} else if (dir.pell_n == -3) {
max_d_n3 = std::max(max_d_n3, dir.d_max);
max_k_n3 = std::max(max_k_n3, dir.k_max);
}
}
const std::vector<std::pair<i64, i64>> pairs_n39 = generate_pell_pairs(-39, max_d_n39, max_k_n39);
const std::vector<std::pair<i64, i64>> pairs_n3 = generate_pell_pairs(-3, max_d_n3, max_k_n3);
const unsigned threads =
choose_thread_count(allow_multithreading, requested_threads, directions.size());
std::vector<std::unordered_map<TriangleKey, u64, TriangleKeyHash>> partial_maps(
static_cast<std::size_t>(threads));
std::vector<std::thread> pool;
pool.reserve(static_cast<std::size_t>(threads));
for (unsigned t = 0U; t < threads; ++t) {
pool.emplace_back([&, t]() {
auto& local = partial_maps[static_cast<std::size_t>(t)];
for (std::size_t i = static_cast<std::size_t>(t); i < directions.size();
i += static_cast<std::size_t>(threads)) {
const Direction& dir = directions[i];
const auto& pairs = (dir.pell_n == -39) ? pairs_n39 : pairs_n3;
for (const auto [d, k_abs] : pairs) {
if (d > dir.d_max || k_abs > dir.k_max) {
continue;
}
for (const i64 sign : {-1LL, 1LL}) {
TriangleKey key;
u64 area = 0ULL;
if (!build_triangle_key(dir, d, sign * k_abs, n, key, area)) {
continue;
}
const auto it = local.find(key);
if (it == local.end()) {
local.emplace(key, area);
} else if (it->second != area) {
std::cerr << "Internal inconsistency: same triangle got two areas.\n";
std::abort();
}
}
}
}
});
}
for (std::thread& th : pool) {
th.join();
}
std::unordered_map<TriangleKey, u64, TriangleKeyHash> merged;
for (auto& local : partial_maps) {
for (auto& kv : local) {
const auto it = merged.find(kv.first);
if (it == merged.end()) {
merged.emplace(std::move(kv));
} else if (it->second != kv.second) {
std::cerr << "Internal inconsistency during merge.\n";
std::abort();
}
}
}
if (out_triangle_count != nullptr) {
*out_triangle_count = merged.size();
}
return sum_triangle_areas_from_map(merged);
}
u128 solve_bruteforce_small(const int n) {
std::unordered_map<TriangleKey, u64, TriangleKeyHash> triangles;
for (int x1 = -n; x1 <= n; ++x1) {
for (int y1 = -n; y1 <= n; ++y1) {
for (int x2 = -n; x2 <= n; ++x2) {
const i64 coeff = static_cast<i64>(x1) + 2LL * static_cast<i64>(x2);
const i64 rhs = -(2LL * static_cast<i64>(x1) + static_cast<i64>(x2)) *
static_cast<i64>(y1);
const auto check_candidate = [&](const i64 y2) {
if (y2 < -n || y2 > n) {
return;
}
const i64 re = static_cast<i64>(x1) * static_cast<i64>(x1) -
static_cast<i64>(y1) * static_cast<i64>(y1) +
static_cast<i64>(x1) * static_cast<i64>(x2) -
static_cast<i64>(y1) * y2 +
static_cast<i64>(x2) * static_cast<i64>(x2) - y2 * y2;
if (re != 39LL) {
return;
}
const i64 x3 = -static_cast<i64>(x1) - static_cast<i64>(x2);
const i64 y3 = -static_cast<i64>(y1) - y2;
if (abs_i64(x3) > n || abs_i64(y3) > n) {
return;
}
std::array<Point, 3> points = {{{x1, y1}, {x2, y2}, {x3, y3}}};
std::sort(points.begin(), points.end());
if (points[0] == points[1] || points[1] == points[2]) {
return;
}
const i64 area2 = abs_i64((points[1].x - points[0].x) * (points[2].y - points[0].y) -
(points[1].y - points[0].y) * (points[2].x - points[0].x));
if (area2 == 0 || (area2 & 1LL) != 0LL) {
return;
}
TriangleKey key;
key.vertices = points;
triangles[key] = static_cast<u64>(area2 / 2LL);
};
if (coeff == 0LL) {
if (rhs != 0LL) {
continue;
}
for (int y2 = -n; y2 <= n; ++y2) {
check_candidate(static_cast<i64>(y2));
}
} else {
if (rhs % coeff != 0LL) {
continue;
}
check_candidate(rhs / coeff);
}
}
}
}
return sum_triangle_areas_from_map(triangles);
}
bool run_checkpoints(const bool allow_multithreading, const unsigned requested_threads) {
struct Checkpoint {
int n;
u64 expected;
};
const std::vector<Checkpoint> known = {
{8, 72ULL},
{10, 252ULL},
{100, 34'632ULL},
{1'000, 3'529'008ULL},
};
for (const Checkpoint& cp : known) {
const u128 got = solve_fast(cp.n, allow_multithreading, requested_threads);
if (got != static_cast<u128>(cp.expected)) {
std::cerr << "Checkpoint failed for A(" << cp.n << "): expected " << cp.expected
<< ", got " << to_string_u128(got) << '\n';
return false;
}
}
const u128 brute20 = solve_bruteforce_small(20);
const u128 fast20 = solve_fast(20, false, 1U);
if (brute20 != fast20) {
std::cerr << "Bruteforce cross-check failed for n=20: brute=" << to_string_u128(brute20)
<< ", fast=" << to_string_u128(fast20) << '\n';
return false;
}
return true;
}
} // namespace
int main(int argc, char** argv) {
std::ios::sync_with_stdio(false);
std::cin.tie(nullptr);
Options options;
if (!parse_arguments(argc, argv, options)) {
return 1;
}
if (options.n < 0) {
std::cerr << "--n must be non-negative.\n";
return 1;
}
if (options.run_checkpoints &&
!run_checkpoints(options.allow_multithreading, options.requested_threads)) {
return 1;
}
std::size_t triangle_count = 0ULL;
const u128 answer = solve_fast(options.n,
options.allow_multithreading,
options.requested_threads,
&triangle_count);
std::cout << to_string_u128(answer) << '\n';
std::cerr << "Triangles counted: " << triangle_count << '\n';
return 0;
}
Python
import math
def generate_pell_pairs(pell_n, d_max, k_max):
if pell_n == -39:
seeds = [(3, 4), (6, 5)]
elif pell_n == -3:
seeds = [(3, 2)]
else:
return []
all_pairs = []
for k0, d0 in seeds:
k, d = k0, d0
while d <= d_max and k <= k_max:
all_pairs.append((d, k))
next_k = 2 * k + 3 * d
next_d = k + 2 * d
k, d = next_k, next_d
return sorted(list(set(all_pairs)))
def build_directions(n):
directions = []
for q0 in (3, 39):
pell_n = -117 // q0
limit = int(math.sqrt(q0)) + 2
for u in range(-limit, limit + 1):
for v in range(-limit, limit + 1):
if u == 0 and v == 0: continue
if math.gcd(abs(u), abs(v)) != 1: continue
if u * u + u * v + v * v != q0: continue
a = 2 * u + v
b = u + 2 * v
m = math.gcd(abs(a), abs(b))
if m != 3: continue
x_scale = max(abs(u), abs(v), abs(u + v))
y_scale_raw = max(abs(a), abs(b), abs(u - v))
if x_scale == 0 or y_scale_raw % m != 0: continue
y_scale = y_scale_raw // m
if y_scale == 0: continue
d_max = n // x_scale
k_max = n // y_scale
if d_max > 0 and k_max > 0:
directions.append({
'q0': q0, 'u': u, 'v': v, 'm': m,
'pell_n': pell_n, 'd_max': d_max, 'k_max': k_max
})
directions.sort(key=lambda d: (d['q0'], d['u'], d['v']))
return directions
def build_triangle_key(dir_obj, d, k, n):
u, v, m = dir_obj['u'], dir_obj['v'], dir_obj['m']
a = 2 * u + v
b = u + 2 * v
if a % m != 0 or b % m != 0 or (u - v) % m != 0:
return None
x1 = d * u
x2 = d * v
x3 = -d * (u + v)
y1 = (b // m) * k
y2 = -(a // m) * k
y3 = ((u - v) // m) * k
if max(abs(x1), abs(x2), abs(x3), abs(y1), abs(y2), abs(y3)) > n:
return None
points = sorted([(x1, y1), (x2, y2), (x3, y3)])
if points[0] == points[1] or points[1] == points[2]:
return None
area = dir_obj['q0'] * abs(d * k)
return tuple(points), area
def solve():
n = 1000000000
directions = build_directions(n)
if not directions: return "0"
max_d_n39 = max_k_n39 = max_d_n3 = max_k_n3 = 0
for d in directions:
if d['pell_n'] == -39:
max_d_n39 = max(max_d_n39, d['d_max'])
max_k_n39 = max(max_k_n39, d['k_max'])
elif d['pell_n'] == -3:
max_d_n3 = max(max_d_n3, d['d_max'])
max_k_n3 = max(max_k_n3, d['k_max'])
pairs_n39 = generate_pell_pairs(-39, max_d_n39, max_k_n39)
pairs_n3 = generate_pell_pairs(-3, max_d_n3, max_k_n3)
triangles = {}
for d_obj in directions:
pairs = pairs_n39 if d_obj['pell_n'] == -39 else pairs_n3
for d, k_abs in pairs:
if d > d_obj['d_max'] or k_abs > d_obj['k_max']:
continue
for sign in (-1, 1):
res = build_triangle_key(d_obj, d, sign * k_abs, n)
if not res: continue
key, area = res
if key not in triangles:
triangles[key] = area
elif triangles[key] != area:
raise Exception("Inconsistency: same triangle got two areas.")
total_area = sum(triangles.values())
return str(total_area)
if __name__ == '__main__':
print(solve())
Java
import java.util.*;
public class Euler385 {
static class Point implements Comparable<Point> {
long x, y;
Point(long x, long y) {
this.x = x;
this.y = y;
}
@Override
public int compareTo(Point o) {
if (this.x != o.x)
return Long.compare(this.x, o.x);
return Long.compare(this.y, o.y);
}
@Override
public boolean equals(Object o) {
if (!(o instanceof Point))
return false;
Point p = (Point) o;
return this.x == p.x && this.y == p.y;
}
@Override
public int hashCode() {
int hx = Long.hashCode(x), hy = Long.hashCode(y);
return hx ^ (hy + 0x9e3779b9 + (hx << 6) + (hx >> 2));
}
}
static class TriangleKey {
Point[] pts;
TriangleKey(Point p1, Point p2, Point p3) {
pts = new Point[] { p1, p2, p3 };
Arrays.sort(pts);
}
@Override
public boolean equals(Object o) {
if (!(o instanceof TriangleKey))
return false;
TriangleKey tk = (TriangleKey) o;
return pts[0].equals(tk.pts[0]) && pts[1].equals(tk.pts[1]) && pts[2].equals(tk.pts[2]);
}
@Override
public int hashCode() {
int h = 0;
for (Point p : pts) {
h ^= p.hashCode() + 0x9e3779b9 + (h << 6) + (h >> 2);
}
return h;
}
}
static class Pair implements Comparable<Pair> {
long d, k;
Pair(long d, long k) {
this.d = d;
this.k = k;
}
@Override
public int compareTo(Pair o) {
if (this.d != o.d)
return Long.compare(this.d, o.d);
return Long.compare(this.k, o.k);
}
@Override
public boolean equals(Object o) {
if (!(o instanceof Pair))
return false;
Pair p = (Pair) o;
return this.d == p.d && this.k == p.k;
}
@Override
public int hashCode() {
return Long.hashCode(d) ^ Long.hashCode(k);
}
}
static class Direction {
int q0, u, v, m, pell_n;
long d_max, k_max;
}
static long gcd(long a, long b) {
while (b != 0) {
long t = a % b;
a = b;
b = t;
}
return Math.abs(a);
}
static List<Pair> generatePellPairs(int pellN, long dMax, long kMax) {
List<Pair> seeds = new ArrayList<>();
if (pellN == -39) {
seeds.add(new Pair(4, 3)); // (d, k) form => k0=3, d0=4
seeds.add(new Pair(5, 6)); // k0=6, d0=5
} else if (pellN == -3) {
seeds.add(new Pair(2, 3)); // k0=3, d0=2
} else {
return new ArrayList<>();
}
List<Pair> all = new ArrayList<>();
for (Pair s : seeds) {
long k = s.k, d = s.d;
while (d <= dMax && k <= kMax) {
all.add(new Pair(d, k));
if (k > (Long.MAX_VALUE - 3 * d) / 2 || d > (Long.MAX_VALUE - k) / 2)
break;
long nK = 2 * k + 3 * d;
long nD = k + 2 * d;
k = nK;
d = nD;
}
}
Set<Pair> set = new HashSet<>(all);
all = new ArrayList<>(set);
Collections.sort(all);
return all;
}
static List<Direction> buildDirections(long n) {
List<Direction> dirs = new ArrayList<>();
int[] q0s = { 3, 39 };
for (int q0 : q0s) {
int pellN = -117 / q0;
int limit = (int) Math.sqrt(q0) + 2;
for (int u = -limit; u <= limit; u++) {
for (int v = -limit; v <= limit; v++) {
if (u == 0 && v == 0)
continue;
if (gcd(u, v) != 1)
continue;
if (u * u + u * v + v * v != q0)
continue;
int a = 2 * u + v;
int b = u + 2 * v;
int m = (int) gcd(a, b);
if (m != 3)
continue;
long xScale = Math.max(Math.max(Math.abs(u), Math.abs(v)), Math.abs(u + v));
long yScaleRaw = Math.max(Math.max(Math.abs(a), Math.abs(b)), Math.abs(u - v));
if (xScale == 0 || yScaleRaw % m != 0)
continue;
long yScale = yScaleRaw / m;
if (yScale == 0)
continue;
long dMax = n / xScale;
long kMax = n / yScale;
if (dMax > 0 && kMax > 0) {
Direction d = new Direction();
d.q0 = q0;
d.u = u;
d.v = v;
d.m = m;
d.pell_n = pellN;
d.d_max = dMax;
d.k_max = kMax;
dirs.add(d);
}
}
}
}
return dirs;
}
static class Result {
TriangleKey tk;
long area;
}
static Result buildTriangleKey(Direction dir, long d, long k, long n) {
long u = dir.u, v = dir.v, m = dir.m;
long a = 2 * u + v;
long b = u + 2 * v;
if (a % m != 0 || b % m != 0 || (u - v) % m != 0)
return null;
long x1 = d * u;
long x2 = d * v;
long x3 = -d * (u + v);
long y1 = (b / m) * k;
long y2 = -(a / m) * k;
long y3 = ((u - v) / m) * k;
if (Math.abs(x1) > n || Math.abs(x2) > n || Math.abs(x3) > n ||
Math.abs(y1) > n || Math.abs(y2) > n || Math.abs(y3) > n) {
return null;
}
Point p1 = new Point(x1, y1), p2 = new Point(x2, y2), p3 = new Point(x3, y3);
TriangleKey tk = new TriangleKey(p1, p2, p3);
if (tk.pts[0].equals(tk.pts[1]) || tk.pts[1].equals(tk.pts[2]))
return null;
Result r = new Result();
r.tk = tk;
r.area = dir.q0 * Math.abs(d * k);
return r;
}
static String solve() {
long n = 1000000000L;
List<Direction> dirs = buildDirections(n);
if (dirs.isEmpty())
return "0";
long maxDN39 = 0, maxKN39 = 0, maxDN3 = 0, maxKN3 = 0;
for (Direction d : dirs) {
if (d.pell_n == -39) {
maxDN39 = Math.max(maxDN39, d.d_max);
maxKN39 = Math.max(maxKN39, d.k_max);
} else if (d.pell_n == -3) {
maxDN3 = Math.max(maxDN3, d.d_max);
maxKN3 = Math.max(maxKN3, d.k_max);
}
}
List<Pair> pairsN39 = generatePellPairs(-39, maxDN39, maxKN39);
List<Pair> pairsN3 = generatePellPairs(-3, maxDN3, maxKN3);
Map<TriangleKey, Long> map = new HashMap<>();
for (Direction dir : dirs) {
List<Pair> pairs = (dir.pell_n == -39) ? pairsN39 : pairsN3;
for (Pair p : pairs) {
if (p.d > dir.d_max || p.k > dir.k_max)
continue;
for (long sign : new long[] { -1, 1 }) {
Result res = buildTriangleKey(dir, p.d, sign * p.k, n);
if (res == null)
continue;
if (map.containsKey(res.tk) && map.get(res.tk) != res.area) {
throw new RuntimeException("Inconsistency in triangle area");
}
map.put(res.tk, res.area);
}
}
}
long total = 0;
for (long a : map.values())
total += a;
return Long.toString(total);
}
public static void main(String[] args) {
System.out.println(solve());
}
}