Problem 397: Triangle on Parabola
View on Project EulerProject Euler Problem 397 Solution
EulerSolve provides an optimized solution for Project Euler Problem 397, Triangle on Parabola, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary For each integer \(k \ge 1\), consider the parabola \(y=\dfrac{x^2}{k}\). We choose three distinct integers \(a \lt b \lt c\) with \(|a|,|b|,|c|\le X\), form the points \(A=\left(a,\dfrac{a^2}{k}\right)\), \(B=\left(b,\dfrac{b^2}{k}\right)\), \(C=\left(c,\dfrac{c^2}{k}\right)\), and ask whether triangle \(ABC\) has at least one \(45^\circ\) angle. Let \(T_k(X)\) be that count for fixed \(k\). The program computes $$F(K,X)=\sum_{k=1}^{K} T_k(X).$$ The key reduction is that chord slopes on this parabola depend only on pair sums such as \(a+b\), so the geometry becomes a divisor problem for \(2k^2\) plus interval counting for the remaining coordinate. Mathematical Approach Step 1: Replace Geometry by Pair Sums Define $$x=a+b,\qquad y=a+c,\qquad z=b+c,$$ so \(x \lt y \lt z\) because \(a \lt b \lt c\). The slope of a chord between two points on \(y=\dfrac{x^2}{k}\) is $$m_{uv}=\frac{v^2-u^2}{k(v-u)}=\frac{u+v}{k}.$$ Therefore the three side slopes are exactly $$m_{AB}=\frac{x}{k},\qquad m_{AC}=\frac{y}{k},\qquad m_{BC}=\frac{z}{k}.$$ This is the crucial simplification: once \((x,y,z)\) are known, the angle conditions no longer involve quadratic expressions in \(a,b,c\). Step 2: Convert the \(45^\circ\) Condition into Diophantine Equations At vertex \(A\), the two scaled direction vectors are proportional to \((k,x)\) and \((k,y)\)....
Detailed mathematical approach
Problem Summary
For each integer \(k \ge 1\), consider the parabola \(y=\dfrac{x^2}{k}\). We choose three distinct integers \(a \lt b \lt c\) with \(|a|,|b|,|c|\le X\), form the points \(A=\left(a,\dfrac{a^2}{k}\right)\), \(B=\left(b,\dfrac{b^2}{k}\right)\), \(C=\left(c,\dfrac{c^2}{k}\right)\), and ask whether triangle \(ABC\) has at least one \(45^\circ\) angle. Let \(T_k(X)\) be that count for fixed \(k\). The program computes
$$F(K,X)=\sum_{k=1}^{K} T_k(X).$$
The key reduction is that chord slopes on this parabola depend only on pair sums such as \(a+b\), so the geometry becomes a divisor problem for \(2k^2\) plus interval counting for the remaining coordinate.
Mathematical Approach
Step 1: Replace Geometry by Pair Sums
Define
$$x=a+b,\qquad y=a+c,\qquad z=b+c,$$
so \(x \lt y \lt z\) because \(a \lt b \lt c\). The slope of a chord between two points on \(y=\dfrac{x^2}{k}\) is
$$m_{uv}=\frac{v^2-u^2}{k(v-u)}=\frac{u+v}{k}.$$
Therefore the three side slopes are exactly
$$m_{AB}=\frac{x}{k},\qquad m_{AC}=\frac{y}{k},\qquad m_{BC}=\frac{z}{k}.$$
This is the crucial simplification: once \((x,y,z)\) are known, the angle conditions no longer involve quadratic expressions in \(a,b,c\).
Step 2: Convert the \(45^\circ\) Condition into Diophantine Equations
At vertex \(A\), the two scaled direction vectors are proportional to \((k,x)\) and \((k,y)\). A \(45^\circ\) angle means that the cross-product magnitude equals the dot product:
$$k(y-x)=k^2+xy.$$
After rearranging, we obtain
$$xy+kx-ky+k^2=0,$$
hence
$$\boxed{(x-k)(y+k)=-2k^2.}$$
At vertex \(B\), the relevant scaled direction vectors are proportional to \((-k,-x)\) and \((k,z)\). The same \(45^\circ\) criterion gives
$$k(z-x)=-(k^2+xz),$$
so
$$\boxed{(x+k)(z-k)=-2k^2.}$$
At vertex \(C\), the symmetric computation yields
$$\boxed{(y-k)(z+k)=-2k^2.}$$
Thus every admissible triangle is encoded by one of three factor equations whose right-hand side is the fixed integer \(-2k^2\). No floating-point trigonometry is needed anywhere in the solver.
Step 3: Enumerate Candidate Pair Sums from Divisors of \(2k^2\)
Write \(pq=-2k^2\). The implementation enumerates every positive divisor \(d\mid 2k^2\), then uses the two signs \(p=\pm d\) and sets \(q=-2k^2/p\).
For the \(A\)- and \(C\)-families we solve \((u-k)(v+k)=-2k^2\), so
$$u=p+k,\qquad v=q-k.$$
For the \(B\)-family we solve \((u+k)(v-k)=-2k^2\), so
$$u=p-k,\qquad v=q+k.$$
Only pairs with \(u \lt v\) and \(|u|,|v|\le 2X\) can come from actual sums of numbers in \([-X,X]\), so the code filters to that range immediately, sorts the resulting lists, and removes duplicates. Different divisors and sign choices can generate the same \((u,v)\), so this deduplication is essential.
Step 4: Count the Missing Coordinate by Integer Intervals
Once a valid pair of sums is known, the third step is not another search; it is just an interval-length computation.
For an angle at \(A\), the pair is \((x,y)=(a+b,a+c)\). Then \(b=x-a\) and \(c=y-a\). Using \(-X\le a,b,c\le X\) and \(a \lt b \lt c\), the allowed values of \(a\) satisfy
$$a\in\left[\max(-X,y-X),\ \min\!\left(X,x+X,\left\lfloor\frac{x-1}{2}\right\rfloor\right)\right].$$
For an angle at \(B\), the pair is \((x,z)=(a+b,b+c)\), so \(a=x-b\), \(c=z-b\), and
$$b\in\left[\max\!\left(-X,z-X,\left\lfloor\frac{x}{2}\right\rfloor+1\right),\ \min\!\left(X,x+X,\left\lfloor\frac{z-1}{2}\right\rfloor\right)\right].$$
For an angle at \(C\), the pair is \((y,z)=(a+c,b+c)\), so \(a=y-c\), \(b=z-c\), and
$$c\in\left[\max\!\left(-X,z-X,\left\lfloor\frac{z}{2}\right\rfloor+1\right),\ \min(X,y+X)\right].$$
The helper functions count_angle_a_for_pair, count_angle_b_for_pair, and
count_angle_c_for_pair implement exactly these formulas. Because sums can be negative, the C++ and
Java versions use a custom floor_div instead of truncating division.
Step 5: Remove Double Counts by Reconstructing \((a,b,c)\)
A triangle can have two \(45^\circ\) angles, so inclusion-exclusion is needed:
$$T_k(X)=N_A+N_B+N_C-N_{AB}-N_{AC}-N_{BC}.$$
When two families agree on the shared sum, we recover the original integer parameters from the ordered sums \(x \lt y \lt z\):
$$a=\frac{x+y-z}{2},\qquad b=\frac{x+z-y}{2},\qquad c=\frac{y+z-x}{2}.$$
The function valid_triangle_from_sums checks exactly what this reconstruction requires: strict order,
even parity of \(x+y+z\), and \(|a|,|b|,|c|\le X\). The overlaps are found by linear merges of sorted pair lists,
so the correction step is \(O(n)\) after sorting.
Worked Example
Take \(k=1\) and \((a,b,c)=(-3,0,2)\). Then the points are \(A=(-3,9)\), \(B=(0,0)\), \(C=(2,4)\), and
$$x=a+b=-3,\qquad z=b+c=2.$$
The middle-angle equation becomes
$$\left(x+1\right)\left(z-1\right)=(-2)(1)=-2=-2\cdot 1^2,$$
so the code classifies this triangle in the \(B\)-family and counts it as a \(45^\circ\) triangle. At a larger checkpoint, the program verifies that
$$F(1,10)=41,$$
which is the same reference value hardcoded in the C++ solution.
How the Code Works
The C++ implementation builds an SPF table once up to the maximum required \(k\), factors each \(2k^2\) in
factor_2k_squared, recursively generates all divisors, constructs the two candidate pair lists in
build_pairs_for_k, evaluates the interval formulas, and subtracts overlaps with three sorted
merge-scans. The C++ version can split the \(k\)-range across threads. The Python file is intentionally not a
second mathematical implementation: it is a thin bridge that compiles and runs the C++ solver, then parses the
final line. The Java file is a direct arithmetic port of the same method and parallelizes with
LongStream.rangeClosed(...).parallel().
Complexity Analysis
Building the SPF sieve up to \(K\) costs \(O(K\log\log K)\) time and \(O(K)\) memory. For fixed \(k\), let \(D_k=d(2k^2)\) be the number of divisors of \(2k^2\). Divisor generation is \(O(D_k)\), interval accumulation is linear in the number of surviving pairs, and sorting/deduplication dominates at roughly \(O(D_k\log D_k)\). The three overlap corrections are linear merges after sorting. Therefore the overall complexity is
$$O\!\left(K\log\log K+\sum_{k=1}^{K} D_k\log D_k\right),$$
with \(O(K)\) global memory for the sieve and \(O(D_k)\) temporary memory per worker.
Footnotes and References
- Problem page: https://projecteuler.net/problem=397
- Parabola and chord geometry: Wikipedia — Parabola
- Angle between two lines: Wikipedia — Angle between two lines
- Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle
- Sieve of Eratosthenes / smallest prime factors: Wikipedia — Sieve of Eratosthenes
Problem 397 source code
C++
#include <algorithm>
#include <cstdint>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <utility>
#include <vector>
namespace {
using i64 = std::int64_t;
using u64 = std::uint64_t;
using u128 = unsigned __int128;
constexpr u64 kDefaultK = 1'000'000ULL;
constexpr i64 kDefaultX = 1'000'000'000LL;
struct Checkpoint {
u64 k = 0ULL;
i64 x = 0LL;
u64 expected = 0ULL;
};
constexpr Checkpoint kCheckpoints[] = {
{1ULL, 10LL, 41ULL},
{10ULL, 100LL, 12'492ULL},
{1'000ULL, 10'000LL, 329'864'366ULL},
};
struct Options {
u64 k = kDefaultK;
i64 x = kDefaultX;
bool run_checkpoints = true;
bool allow_multithreading = true;
unsigned requested_threads = 0U;
};
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_i64_after_prefix(const std::string& arg, const char* prefix, i64& 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;
}
if (parsed > static_cast<u64>(std::numeric_limits<i64>::max())) {
return false;
}
value = static_cast<i64>(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, "--k=", parsed_u64)) {
options.k = parsed_u64;
continue;
}
i64 parsed_i64 = 0LL;
if (parse_i64_after_prefix(arg, "--x=", parsed_i64)) {
options.x = parsed_i64;
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;
}
if (options.k == 0ULL) {
std::cerr << "--k must be positive.\n";
return false;
}
if (options.x <= 0LL) {
std::cerr << "--x must be positive.\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 = 1'000ULL;
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 floor_div(const i64 a, const i64 b) {
if (a >= 0LL) {
return a / b;
}
return -(((-a) + b - 1LL) / b);
}
u64 count_interval(const i64 lo, const i64 hi) {
if (lo > hi) {
return 0ULL;
}
return static_cast<u64>(hi - lo + 1LL);
}
u64 count_angle_a_for_pair(const i64 u, const i64 v, const i64 x_bound) {
const i64 lo = std::max(-x_bound, v - x_bound);
const i64 hi = std::min({x_bound, u + x_bound, floor_div(u - 1LL, 2LL)});
return count_interval(lo, hi);
}
u64 count_angle_b_for_pair(const i64 u, const i64 v, const i64 x_bound) {
const i64 lo = std::max({-x_bound, v - x_bound, floor_div(u, 2LL) + 1LL});
const i64 hi = std::min({x_bound, u + x_bound, floor_div(v - 1LL, 2LL)});
return count_interval(lo, hi);
}
u64 count_angle_c_for_pair(const i64 u, const i64 v, const i64 x_bound) {
const i64 lo = std::max({-x_bound, v - x_bound, floor_div(v, 2LL) + 1LL});
const i64 hi = std::min(x_bound, u + x_bound);
return count_interval(lo, hi);
}
bool valid_triangle_from_sums(const i64 s1, const i64 s2, const i64 s3, const i64 x_bound) {
if (!(s1 < s2 && s2 < s3)) {
return false;
}
const i64 parity_sum = s1 + s2 + s3;
if ((parity_sum & 1LL) != 0LL) {
return false;
}
const i64 a = (s1 + s2 - s3) / 2LL;
const i64 b = (s1 + s3 - s2) / 2LL;
const i64 c = (s2 + s3 - s1) / 2LL;
if (a < -x_bound || a > x_bound) {
return false;
}
if (b < -x_bound || b > x_bound) {
return false;
}
if (c < -x_bound || c > x_bound) {
return false;
}
return a < b && b < c;
}
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<int> build_spf(const u64 limit) {
std::vector<int> spf(static_cast<std::size_t>(limit) + 1ULL, 0);
if (limit >= 1ULL) {
spf[1] = 1;
}
for (u64 i = 2ULL; i <= limit; ++i) {
if (spf[static_cast<std::size_t>(i)] == 0) {
spf[static_cast<std::size_t>(i)] = static_cast<int>(i);
if (i <= limit / i) {
for (u64 j = i * i; j <= limit; j += i) {
if (spf[static_cast<std::size_t>(j)] == 0) {
spf[static_cast<std::size_t>(j)] = static_cast<int>(i);
}
}
}
}
}
return spf;
}
void collect_divisors_recursive(const std::vector<std::pair<u64, int>>& factors,
const std::size_t index,
const u64 current,
std::vector<u64>& divisors) {
if (index == factors.size()) {
divisors.push_back(current);
return;
}
const u64 prime = factors[index].first;
const int exponent = factors[index].second;
u64 value = current;
for (int e = 0; e <= exponent; ++e) {
collect_divisors_recursive(factors, index + 1ULL, value, divisors);
if (e != exponent) {
value *= prime;
}
}
}
void factor_2k_squared(const u64 k,
const std::vector<int>& spf,
std::vector<std::pair<u64, int>>& factors) {
factors.clear();
u64 n = k;
int exponent_two_in_k = 0;
while ((n & 1ULL) == 0ULL) {
++exponent_two_in_k;
n >>= 1U;
}
factors.emplace_back(2ULL, 2 * exponent_two_in_k + 1);
while (n > 1ULL) {
const int p = spf[static_cast<std::size_t>(n)];
int exp_in_k = 0;
while (n % static_cast<u64>(p) == 0ULL) {
n /= static_cast<u64>(p);
++exp_in_k;
}
factors.emplace_back(static_cast<u64>(p), 2 * exp_in_k);
}
}
void build_pairs_for_k(const u64 k,
const i64 x_bound,
const std::vector<u64>& divisors,
std::vector<std::pair<i64, i64>>& plus_pairs,
std::vector<std::pair<i64, i64>>& minus_pairs,
std::vector<std::pair<i64, i64>>& plus_out,
std::vector<std::pair<i64, i64>>& plus_in,
std::vector<std::pair<i64, i64>>& minus_out,
std::vector<std::pair<i64, i64>>& minus_in) {
plus_pairs.clear();
minus_pairs.clear();
plus_out.clear();
plus_in.clear();
minus_out.clear();
minus_in.clear();
const i64 kk = static_cast<i64>(k);
const i64 n = static_cast<i64>(2ULL * k * k);
const i64 sum_bound = 2LL * x_bound;
plus_pairs.reserve(divisors.size() * 2ULL);
minus_pairs.reserve(divisors.size() * 2ULL);
for (const u64 d_u64 : divisors) {
const i64 d = static_cast<i64>(d_u64);
for (const int sign : {-1, 1}) {
const i64 p = static_cast<i64>(sign) * d;
const i64 q = -n / p;
const i64 u_plus = p + kk;
const i64 v_plus = q - kk;
if (u_plus < v_plus &&
-sum_bound <= u_plus && u_plus <= sum_bound &&
-sum_bound <= v_plus && v_plus <= sum_bound) {
plus_pairs.emplace_back(u_plus, v_plus);
}
const i64 u_minus = p - kk;
const i64 v_minus = q + kk;
if (u_minus < v_minus &&
-sum_bound <= u_minus && u_minus <= sum_bound &&
-sum_bound <= v_minus && v_minus <= sum_bound) {
minus_pairs.emplace_back(u_minus, v_minus);
}
}
}
std::sort(plus_pairs.begin(), plus_pairs.end());
plus_pairs.erase(std::unique(plus_pairs.begin(), plus_pairs.end()), plus_pairs.end());
std::sort(minus_pairs.begin(), minus_pairs.end());
minus_pairs.erase(std::unique(minus_pairs.begin(), minus_pairs.end()), minus_pairs.end());
plus_out = plus_pairs;
minus_out = minus_pairs;
plus_in.reserve(plus_pairs.size());
for (const auto& [u, v] : plus_pairs) {
plus_in.emplace_back(v, u);
}
std::sort(plus_in.begin(), plus_in.end());
minus_in.reserve(minus_pairs.size());
for (const auto& [u, v] : minus_pairs) {
minus_in.emplace_back(v, u);
}
std::sort(minus_in.begin(), minus_in.end());
}
u128 solve_range(const u64 thread_index,
const u64 thread_count,
const u64 k_limit,
const i64 x_bound,
const std::vector<int>& spf) {
u128 subtotal = 0U;
std::vector<std::pair<u64, int>> factors;
std::vector<u64> divisors;
std::vector<std::pair<i64, i64>> plus_pairs;
std::vector<std::pair<i64, i64>> minus_pairs;
std::vector<std::pair<i64, i64>> plus_out;
std::vector<std::pair<i64, i64>> plus_in;
std::vector<std::pair<i64, i64>> minus_out;
std::vector<std::pair<i64, i64>> minus_in;
factors.reserve(16ULL);
divisors.reserve(2048ULL);
for (u64 k = 1ULL + thread_index; k <= k_limit; k += thread_count) {
factor_2k_squared(k, spf, factors);
divisors.clear();
collect_divisors_recursive(factors, 0ULL, 1ULL, divisors);
build_pairs_for_k(k,
x_bound,
divisors,
plus_pairs,
minus_pairs,
plus_out,
plus_in,
minus_out,
minus_in);
u128 count_a = 0U;
u128 count_b = 0U;
u128 count_c = 0U;
for (const auto& [u, v] : plus_out) {
count_a += static_cast<u128>(count_angle_a_for_pair(u, v, x_bound));
count_c += static_cast<u128>(count_angle_c_for_pair(u, v, x_bound));
}
for (const auto& [u, v] : minus_out) {
count_b += static_cast<u128>(count_angle_b_for_pair(u, v, x_bound));
}
u128 overlap_ab = 0U;
{
std::size_t i = 0ULL;
std::size_t j = 0ULL;
while (i < plus_out.size() && j < minus_out.size()) {
const i64 key_plus = plus_out[i].first;
const i64 key_minus = minus_out[j].first;
if (key_plus < key_minus) {
++i;
} else if (key_minus < key_plus) {
++j;
} else {
const i64 x = key_plus;
const i64 y = plus_out[i].second;
const i64 z = minus_out[j].second;
if (valid_triangle_from_sums(x, y, z, x_bound)) {
++overlap_ab;
}
++i;
++j;
}
}
}
u128 overlap_ac = 0U;
{
std::size_t i = 0ULL;
std::size_t j = 0ULL;
while (i < plus_in.size() && j < plus_out.size()) {
const i64 key_in = plus_in[i].first;
const i64 key_out = plus_out[j].first;
if (key_in < key_out) {
++i;
} else if (key_out < key_in) {
++j;
} else {
const i64 y = key_in;
const i64 x = plus_in[i].second;
const i64 z = plus_out[j].second;
if (valid_triangle_from_sums(x, y, z, x_bound)) {
++overlap_ac;
}
++i;
++j;
}
}
}
u128 overlap_bc = 0U;
{
std::size_t i = 0ULL;
std::size_t j = 0ULL;
while (i < plus_in.size() && j < minus_in.size()) {
const i64 key_plus = plus_in[i].first;
const i64 key_minus = minus_in[j].first;
if (key_plus < key_minus) {
++i;
} else if (key_minus < key_plus) {
++j;
} else {
const i64 z = key_plus;
const i64 y = plus_in[i].second;
const i64 x = minus_in[j].second;
if (valid_triangle_from_sums(x, y, z, x_bound)) {
++overlap_bc;
}
++i;
++j;
}
}
}
subtotal += count_a + count_b + count_c - overlap_ab - overlap_ac - overlap_bc;
}
return subtotal;
}
u128 solve_f(const u64 k_limit,
const i64 x_bound,
const bool allow_multithreading,
const unsigned requested_threads,
const std::vector<int>& spf) {
const unsigned thread_count =
choose_thread_count(allow_multithreading, requested_threads, static_cast<std::size_t>(k_limit));
if (thread_count == 1U) {
return solve_range(0ULL, 1ULL, k_limit, x_bound, spf);
}
std::vector<u128> partial(thread_count, 0U);
std::vector<std::thread> workers;
workers.reserve(thread_count);
for (unsigned t = 0U; t < thread_count; ++t) {
workers.emplace_back([&, t]() {
partial[t] = solve_range(static_cast<u64>(t),
static_cast<u64>(thread_count),
k_limit,
x_bound,
spf);
});
}
for (std::thread& worker : workers) {
worker.join();
}
u128 total = 0U;
for (const u128 value : partial) {
total += value;
}
return total;
}
} // namespace
int main(int argc, char** argv) {
Options options;
if (!parse_arguments(argc, argv, options)) {
std::cerr << "Usage: " << argv[0]
<< " [--k=<value>] [--x=<value>] [--threads=<value>]"
<< " [--single-thread] [--skip-checkpoints]\n";
return 1;
}
u64 max_k_needed = options.k;
if (options.run_checkpoints) {
for (const Checkpoint& cp : kCheckpoints) {
max_k_needed = std::max(max_k_needed, cp.k);
}
}
const std::vector<int> spf = build_spf(max_k_needed);
if (options.run_checkpoints) {
for (const Checkpoint& cp : kCheckpoints) {
const u128 actual =
solve_f(cp.k, cp.x, options.allow_multithreading, options.requested_threads, spf);
const u128 expected = static_cast<u128>(cp.expected);
std::cout << "Checkpoint F(" << cp.k << ", " << cp.x << ") = "
<< to_string_u128(actual) << " (expected " << cp.expected << ")\n";
if (actual != expected) {
std::cerr << "Checkpoint mismatch.\n";
return 1;
}
}
}
const u128 answer =
solve_f(options.k, options.x, options.allow_multithreading, options.requested_threads, spf);
std::cout << "F(" << options.k << ", " << options.x << ") = " << to_string_u128(answer)
<< '\n';
return 0;
}
Python
from __future__ import annotations
import re
import shutil
import subprocess
from pathlib import Path
ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")
def parse_output(stdout: str) -> str:
lines = [line.strip() for line in stdout.splitlines() if line.strip()]
if not lines:
return ""
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.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.List;
import java.util.stream.LongStream;
public class Euler397 {
static long floorDiv(long a, long b) {
if (a >= 0)
return a / b;
return -(((-a) + b - 1) / b);
}
static long countInterval(long lo, long hi) {
if (lo > hi)
return 0;
return hi - lo + 1;
}
static long countAngleA(long u, long v, long xBound) {
long lo = Math.max(-xBound, v - xBound);
long hi = Math.min(xBound, Math.min(u + xBound, floorDiv(u - 1, 2)));
return countInterval(lo, hi);
}
static long countAngleB(long u, long v, long xBound) {
long lo = Math.max(-xBound, Math.max(v - xBound, floorDiv(u, 2) + 1));
long hi = Math.min(xBound, Math.min(u + xBound, floorDiv(v - 1, 2)));
return countInterval(lo, hi);
}
static long countAngleC(long u, long v, long xBound) {
long lo = Math.max(-xBound, Math.max(v - xBound, floorDiv(v, 2) + 1));
long hi = Math.min(xBound, u + xBound);
return countInterval(lo, hi);
}
static boolean validTriangle(long s1, long s2, long s3, long xBound) {
if (!(s1 < s2 && s2 < s3))
return false;
long parity = s1 + s2 + s3;
if ((parity & 1) != 0)
return false;
long a = (s1 + s2 - s3) / 2;
long b = (s1 + s3 - s2) / 2;
long c = (s2 + s3 - s1) / 2;
if (a < -xBound || a > xBound)
return false;
if (b < -xBound || b > xBound)
return false;
if (c < -xBound || c > xBound)
return false;
return a < b && b < c;
}
static class Factor {
long p;
int e;
Factor(long p, int e) {
this.p = p;
this.e = e;
}
}
static class Pair implements Comparable<Pair> {
long first, second;
Pair(long f, long s) {
first = f;
second = s;
}
@Override
public int compareTo(Pair o) {
if (this.first != o.first)
return Long.compare(this.first, o.first);
return Long.compare(this.second, o.second);
}
@Override
public boolean equals(Object o) {
if (!(o instanceof Pair))
return false;
Pair p = (Pair) o;
return this.first == p.first && this.second == p.second;
}
@Override
public int hashCode() {
return Long.hashCode(first) ^ Long.hashCode(second);
}
}
static void getDivisors(List<Factor> factors, int idx, long current, List<Long> result) {
if (idx == factors.size()) {
result.add(current);
return;
}
long p = factors.get(idx).p;
int exp = factors.get(idx).e;
long val = current;
for (int i = 0; i <= exp; i++) {
getDivisors(factors, idx + 1, val, result);
if (i != exp)
val *= p;
}
}
static long solveForK(long k, long xBound, int[] spf) {
long n_k = k;
int exp2 = 0;
while ((n_k & 1) == 0) {
exp2++;
n_k >>= 1;
}
List<Factor> factors = new ArrayList<>();
factors.add(new Factor(2, 2 * exp2 + 1));
while (n_k > 1) {
int p = spf[(int) n_k];
int expP = 0;
while (n_k % p == 0) {
n_k /= p;
expP++;
}
factors.add(new Factor(p, 2 * expP));
}
List<Long> divisors = new ArrayList<>();
getDivisors(factors, 0, 1L, divisors);
long kk = k;
long nVal = 2 * k * k;
long sumBound = 2 * xBound;
List<Pair> plusOut = new ArrayList<>();
List<Pair> minusOut = new ArrayList<>();
for (long d : divisors) {
for (int sign : new int[] { -1, 1 }) {
long p = sign * d;
long q = -nVal / p;
long uPlus = p + kk;
long vPlus = q - kk;
if (uPlus < vPlus && uPlus >= -sumBound && uPlus <= sumBound && vPlus >= -sumBound
&& vPlus <= sumBound) {
plusOut.add(new Pair(uPlus, vPlus));
}
long uMinus = p - kk;
long vMinus = q + kk;
if (uMinus < vMinus && uMinus >= -sumBound && uMinus <= sumBound && vMinus >= -sumBound
&& vMinus <= sumBound) {
minusOut.add(new Pair(uMinus, vMinus));
}
}
}
Collections.sort(plusOut);
List<Pair> plusOutUnique = new ArrayList<>();
if (!plusOut.isEmpty()) {
plusOutUnique.add(plusOut.get(0));
for (int i = 1; i < plusOut.size(); i++) {
if (!plusOut.get(i).equals(plusOutUnique.get(plusOutUnique.size() - 1))) {
plusOutUnique.add(plusOut.get(i));
}
}
}
plusOut = plusOutUnique;
Collections.sort(minusOut);
List<Pair> minusOutUnique = new ArrayList<>();
if (!minusOut.isEmpty()) {
minusOutUnique.add(minusOut.get(0));
for (int i = 1; i < minusOut.size(); i++) {
if (!minusOut.get(i).equals(minusOutUnique.get(minusOutUnique.size() - 1))) {
minusOutUnique.add(minusOut.get(i));
}
}
}
minusOut = minusOutUnique;
List<Pair> plusIn = new ArrayList<>(plusOut.size());
for (Pair p : plusOut)
plusIn.add(new Pair(p.second, p.first));
Collections.sort(plusIn);
List<Pair> minusIn = new ArrayList<>(minusOut.size());
for (Pair p : minusOut)
minusIn.add(new Pair(p.second, p.first));
Collections.sort(minusIn);
long countA = 0, countB = 0, countC = 0;
for (Pair p : plusOut) {
countA += countAngleA(p.first, p.second, xBound);
countC += countAngleC(p.first, p.second, xBound);
}
for (Pair p : minusOut) {
countB += countAngleB(p.first, p.second, xBound);
}
long overlapAB = 0;
int i = 0, j = 0;
while (i < plusOut.size() && j < minusOut.size()) {
long keyP = plusOut.get(i).first;
long keyM = minusOut.get(j).first;
if (keyP < keyM)
i++;
else if (keyM < keyP)
j++;
else {
long x = keyP;
long y = plusOut.get(i).second;
long z = minusOut.get(j).second;
if (validTriangle(x, y, z, xBound))
overlapAB++;
i++;
j++;
}
}
long overlapAC = 0;
i = 0;
j = 0;
while (i < plusIn.size() && j < plusOut.size()) {
long keyIn = plusIn.get(i).first;
long keyOut = plusOut.get(j).first;
if (keyIn < keyOut)
i++;
else if (keyOut < keyIn)
j++;
else {
long y = keyIn;
long x = plusIn.get(i).second;
long z = plusOut.get(j).second;
if (validTriangle(x, y, z, xBound))
overlapAC++;
i++;
j++;
}
}
long overlapBC = 0;
i = 0;
j = 0;
while (i < plusIn.size() && j < minusIn.size()) {
long keyP = plusIn.get(i).first;
long keyM = minusIn.get(j).first;
if (keyP < keyM)
i++;
else if (keyM < keyP)
j++;
else {
long z = keyP;
long y = plusIn.get(i).second;
long x = minusIn.get(j).second;
if (validTriangle(x, y, z, xBound))
overlapBC++;
i++;
j++;
}
}
return countA + countB + countC - overlapAB - overlapAC - overlapBC;
}
static String solve() {
long kLimit = 1000000L;
long xBound = 1000000000L;
int[] spf = new int[(int) kLimit + 1];
for (int i = 2; i <= kLimit; i++) {
if (spf[i] == 0) {
spf[i] = i;
if ((long) i * i <= kLimit) {
for (int j = i * i; j <= kLimit; j += i) {
if (spf[j] == 0)
spf[j] = i;
}
}
}
}
long total = LongStream.rangeClosed(1, kLimit)
.parallel()
.map(k -> solveForK(k, xBound, spf))
.sum();
return Long.toString(total);
}
public static void main(String[] args) {
System.out.println(solve());
}
}