Problem 579: Lattice Points in Lattice Cubes

Project Euler Problem 579 Solution

EulerSolve provides an optimized solution for Project Euler Problem 579, Lattice Points in Lattice Cubes, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $S(n)$ denote the sum, over every lattice cube that can be translated to lie inside $[0,n]^3$, of the number of lattice points contained in that cube. A lattice cube means that all eight vertices lie in $\mathbb{Z}^3$. A direct search over all vertex sets would be hopelessly redundant. The implementation instead enumerates primitive cube frames, extends them by every admissible integer scale, counts how many integer translations fit inside the box, and adds the exact lattice-point count of each scaled cube. Mathematical Approach Describe one cube by three edge vectors $u,v,w\in\mathbb{Z}^3$ issuing from a single vertex. They must satisfy $$u\cdot v=u\cdot w=v\cdot w=0,\qquad u\cdot u=v\cdot v=w\cdot w.$$ The key is that such integer triples can be generated systematically from integer quaternions, and once one primitive frame is known, all larger cubes in the same orientation follow by integer scaling. Step 1: Generate primitive orthogonal edge triples from integer quaternions Start from integers $a,b,c,d$, not all zero, and define $$\begin{aligned} u&=\frac{1}{g}\left(a^2+b^2-c^2-d^2,\ 2(bc-ad),\ 2(bd+ac)\right),\\ v&=\frac{1}{g}\left(2(bc+ad),\ a^2-b^2+c^2-d^2,\ 2(cd-ab)\right),\\ w&=\frac{1}{g}\left(2(bd-ac),\ 2(cd+ab),\ a^2-b^2-c^2+d^2\right)....

Detailed mathematical approach

Problem Summary

Let $S(n)$ denote the sum, over every lattice cube that can be translated to lie inside $[0,n]^3$, of the number of lattice points contained in that cube. A lattice cube means that all eight vertices lie in $\mathbb{Z}^3$.

A direct search over all vertex sets would be hopelessly redundant. The implementation instead enumerates primitive cube frames, extends them by every admissible integer scale, counts how many integer translations fit inside the box, and adds the exact lattice-point count of each scaled cube.

Mathematical Approach

Describe one cube by three edge vectors $u,v,w\in\mathbb{Z}^3$ issuing from a single vertex. They must satisfy

$$u\cdot v=u\cdot w=v\cdot w=0,\qquad u\cdot u=v\cdot v=w\cdot w.$$

The key is that such integer triples can be generated systematically from integer quaternions, and once one primitive frame is known, all larger cubes in the same orientation follow by integer scaling.

Step 1: Generate primitive orthogonal edge triples from integer quaternions

Start from integers $a,b,c,d$, not all zero, and define

$$\begin{aligned} u&=\frac{1}{g}\left(a^2+b^2-c^2-d^2,\ 2(bc-ad),\ 2(bd+ac)\right),\\ v&=\frac{1}{g}\left(2(bc+ad),\ a^2-b^2+c^2-d^2,\ 2(cd-ab)\right),\\ w&=\frac{1}{g}\left(2(bd-ac),\ 2(cd+ab),\ a^2-b^2-c^2+d^2\right). \end{aligned}$$

The divisor $g$ removes the compulsory common factor forced by parity:

$$g=\begin{cases} 4,& \text{if }a,b,c,d\text{ are all odd},\\ 2,& \text{if exactly two of them are odd},\\ 1,& \text{otherwise}. \end{cases}$$

When $\gcd(a,b,c,d)=1$, these vectors are primitive in the sense needed by the search, and they satisfy

$$u\cdot v=u\cdot w=v\cdot w=0,\qquad u\cdot u=v\cdot v=w\cdot w=\left(\frac{a^2+b^2+c^2+d^2}{g}\right)^2.$$

If we set

$$N_0=a^2+b^2+c^2+d^2,\qquad L_0=\frac{N_0}{g},$$

then $L_0$ is the Euclidean edge length of the primitive cube frame.

Step 2: Pass from one primitive frame to every integer scale

For any integer scale $s\ge 1$, the same orientation gives a larger cube with edges

$$su,\qquad sv,\qquad sw,$$

and edge length

$$L=sL_0.$$

For a lattice vector $x=(x_1,x_2,x_3)$, write

$$\delta(x)=\gcd(|x_1|,|x_2|,|x_3|).$$

The segment from $0$ to $sx$ then contains $s\,\delta(x)+1$ lattice points. For the three primitive edges define

$$\Delta_0=\delta(u)+\delta(v)+\delta(w),\qquad \Delta=s\Delta_0.$$

These quantities measure the lattice step size along the three edge directions and are exactly what the point-counting polynomial needs.

Step 3: Derive the exact lattice-point count of one scaled cube

A standard lattice-parallelepiped formula says that for edge vectors $p,q,r\in\mathbb{Z}^3$,

$$\#\bigl((\text{parallelepiped}(p,q,r))\cap\mathbb{Z}^3\bigr)=|\det(p,q,r)|+\delta(p\times q)+\delta(q\times r)+\delta(r\times p)+\delta(p)+\delta(q)+\delta(r)+1.$$

Apply this with $p=su$, $q=sv$, $r=sw$. Because the edges are mutually orthogonal and all have length $L$, we get

$$|\det(su,sv,sw)|=L^3.$$

Also, orthogonality and equal edge lengths imply

$$ (su)\times(sv)=\pm L(sw),\qquad (sv)\times(sw)=\pm L(su),\qquad (sw)\times(su)=\pm L(sv).$$

Therefore

$$\delta((su)\times(sv))=L\,\delta(sw),\qquad \delta((sv)\times(sw))=L\,\delta(su),\qquad \delta((sw)\times(su))=L\,\delta(sv).$$

Summing the three face terms gives $L\Delta$, and summing the edge terms gives $\Delta$. Hence the lattice-point count of one scaled cube is

$$\boxed{\operatorname{points}(s)=L^3+(L+1)\Delta+1.}$$

As a sanity check, for the axis-aligned cube of side $m$ we have $L=m$ and $\Delta=3m$, so the formula becomes $(m+1)^3$, exactly as expected.

Step 4: Count how many translations fit inside $[0,n]^3$

For the primitive frame define the axis spans

$$\sigma_x=|u_1|+|v_1|+|w_1|,\qquad \sigma_y=|u_2|+|v_2|+|w_2|,\qquad \sigma_z=|u_3|+|v_3|+|w_3|.$$

These are the side lengths of the smallest axis-aligned box containing the primitive cube. The reason is that each vertex coordinate is a subset sum of the edge coordinates, so the total range along one axis is the sum of the absolute values on that axis.

After scaling by $s$, the spans become $s\sigma_x$, $s\sigma_y$, $s\sigma_z$. The number of admissible integer translations is therefore

$$T(s)=(n+1-s\sigma_x)(n+1-s\sigma_y)(n+1-s\sigma_z),$$

provided every factor is positive. Equivalently, the maximum allowed scale is

$$s_{\max}=\min\left(\left\lfloor\frac{n}{\sigma_x}\right\rfloor,\left\lfloor\frac{n}{\sigma_y}\right\rfloor,\left\lfloor\frac{n}{\sigma_z}\right\rfloor\right)=\left\lfloor\frac{n}{\max(\sigma_x,\sigma_y,\sigma_z)}\right\rfloor.$$

Step 5: Sum every frame, every scale, then divide by the 24-fold symmetry

For one primitive frame, the raw contribution is

$$\sum_{s=1}^{s_{\max}} T(s)\left((sL_0)^3+(sL_0+1)s\Delta_0+1\right).$$

The quaternion parameterization enumerates ordered orthogonal edge triples, not unlabeled geometric cubes. The same geometric cube appears in $24$ equivalent orientations, corresponding to the rotational symmetry group of the cube. So if $R(n)$ is the total raw sum over all primitive frames, then

$$\boxed{S(n)=\frac{R(n)}{24}.}$$

This is why the implementation checks that the raw total is divisible by $24$ before reporting the final answer.

Step 6: Worked example for $n=2$

Take the primitive axis frame

$$u=(1,0,0),\qquad v=(0,1,0),\qquad w=(0,0,1).$$

Then $L_0=1$, $\Delta_0=3$, and $\sigma_x=\sigma_y=\sigma_z=1$, so $s_{\max}=2$.

For $s=1$, we have

$$T(1)=2^3=8,\qquad \operatorname{points}(1)=1^3+(1+1)\cdot 3+1=8,$$

giving a contribution of $8\cdot 8=64$.

For $s=2$, we have

$$T(2)=1,\qquad \operatorname{points}(2)=2^3+(2+1)\cdot 6+1=27,$$

giving a contribution of $27$.

No larger scale fits, so

$$S(2)=64+27=91,$$

which matches the checkpoint used by the implementations.

How the Code Works

The C++, Python, and Java implementations all use the same underlying search. The core solver splits the quaternion enumeration by parity pattern, because the parity determines the normalization divisor $g$ and the norm bound $N_0\le gn$. That lets it prebuild only the integer values that can actually occur in each case.

During enumeration, partial sums of $a^2+b^2+c^2+d^2$ are checked early, so branches that already exceed the bound are abandoned immediately. The search also rejects non-primitive quadruples by gcd tests and uses canonical sign restrictions to avoid the trivial duplication coming from replacing the parameter quadruple by its negation.

For every surviving frame, the implementation computes the three edge vectors, rejects frames whose axis spans already exceed the box, derives $L_0$, $\Delta_0$, and $s_{\max}$, and then loops over all valid scales $s$. At each scale it multiplies the translation count by the cube point-count formula, adds the result to a thread-local accumulator, and only after all threads finish divides the raw total by $24$.

The Python and Java versions are thin wrappers around the same compiled solver, so all three languages share the same mathematics, the same checkpoints, and the same final value.

Complexity Analysis

The dominant cost is the quaternion search. A loose geometric upper bound comes from the four-dimensional ball $a^2+b^2+c^2+d^2\le 4n$, which contains $O(n^2)$ integer tuples, so the worst-case search space is roughly quadratic in $n$.

In practice the implementation is much faster than that crude bound suggests. Partial norm checks stop many branches early, gcd filters remove non-primitive tuples, span checks discard frames that can never fit in the box, and each surviving frame only iterates up to $s_{\max}$. The method is therefore strongly prune-driven rather than a naive enumeration of all cubes.

Memory usage is modest. The precomputed even and odd value tables only reach $O(\sqrt{n})$ in magnitude, and the parallel version stores only task metadata plus one accumulator per worker thread, so the overall memory footprint stays small compared with the arithmetic work.

Footnotes and References

Problem page: Project Euler 579
Quaternion background: Wikipedia - Quaternion
Euler-Rodrigues parameterization: Wikipedia - Euler-Rodrigues formula
Lattice polytopes: Wikipedia - Lattice polytope
Ehrhart theory: Wikipedia - Ehrhart polynomial

Problem 579 source code

C++

#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <numeric>
#include <thread>
#include <vector>

using namespace std;

namespace {

struct ValList {
    vector<int> v;
    vector<int> v2;
};

struct ValueSets {
    ValList even;
    ValList odd;
    ValList even_nonneg;
    ValList odd_nonneg;
};

struct Pattern {
    const ValList* a_vals;
    const ValList* b_vals;
    const ValList* c_vals;
    const ValList* d_vals;
    int g;
    int nmax;
};

struct Task {
    int pattern_idx;
    int a_start;
    int a_end;
};

inline int64_t iabs(int64_t x) {
    return x < 0 ? -x : x;
}

inline int gcd3(int64_t a, int64_t b, int64_t c) {
    return std::gcd((int)iabs(a), std::gcd((int)iabs(b), (int)iabs(c)));
}

ValueSets build_values(int M) {
    ValueSets out;
    out.even.v.reserve(M + 1);
    out.odd.v.reserve(M + 1);
    out.even_nonneg.v.reserve(M + 1);
    out.odd_nonneg.v.reserve(M + 1);

    for (int k = 0; k <= M; ++k) {
        if ((k & 1) == 0) {
            if (k == 0) {
                out.even.v.push_back(0);
                out.even_nonneg.v.push_back(0);
            } else {
                out.even.v.push_back(k);
                out.even.v.push_back(-k);
                out.even_nonneg.v.push_back(k);
            }
        } else {
            out.odd.v.push_back(k);
            out.odd.v.push_back(-k);
            out.odd_nonneg.v.push_back(k);
        }
    }

    auto fill_sq = [](ValList& list) {
        list.v2.resize(list.v.size());
        for (size_t i = 0; i < list.v.size(); ++i) {
            int x = list.v[i];
            list.v2[i] = x * x;
        }
    };

    fill_sq(out.even);
    fill_sq(out.odd);
    fill_sq(out.even_nonneg);
    fill_sq(out.odd_nonneg);
    return out;
}

unsigned __int128 compute_sum_frames(int n, unsigned int threads) {
    int m1 = (int)sqrt((long double)n);
    int m2 = (int)sqrt((long double)(2LL * n));
    int m3 = (int)sqrt((long double)(4LL * n));

    ValueSets vals1 = build_values(m1);
    ValueSets vals2 = build_values(m2);
    ValueSets vals3 = build_values(m3);

    vector<Pattern> patterns;
    patterns.reserve(15);

    for (int mask = 1; mask < 16; ++mask) {
        int odd_count = __builtin_popcount((unsigned)mask);
        if (odd_count == 0) continue;

        // For primitive quaternions, the frame gcd is 1/2/4 based on odd-component count.
        int g = 1;
        if (odd_count == 2) g = 2;
        else if (odd_count == 4) g = 4;

        int nmax = n * g;

        ValueSets* vs = nullptr;
        if (g == 1) vs = &vals1;
        else if (g == 2) vs = &vals2;
        else vs = &vals3;

        auto pick_full = [&](int bit) -> const ValList* {
            return (bit ? &vs->odd : &vs->even);
        };
        auto pick_nonneg = [&](int bit) -> const ValList* {
            return (bit ? &vs->odd_nonneg : &vs->even_nonneg);
        };

        Pattern pat;
        pat.a_vals = pick_nonneg(mask & 1);
        pat.b_vals = pick_full(mask & 2);
        pat.c_vals = pick_full(mask & 4);
        pat.d_vals = pick_full(mask & 8);
        pat.g = g;
        pat.nmax = nmax;
        patterns.push_back(pat);
    }

    vector<Task> tasks;
    const int chunk = 4;
    for (int idx = 0; idx < (int)patterns.size(); ++idx) {
        int a_size = (int)patterns[idx].a_vals->v.size();
        for (int i = 0; i < a_size; i += chunk) {
            Task t{idx, i, min(i + chunk, a_size)};
            tasks.push_back(t);
        }
    }

    if (threads == 0) threads = 1;
    vector<unsigned __int128> sums(threads, 0);
    atomic<size_t> task_index{0};

    auto worker = [&](unsigned int tid) {
        unsigned __int128 local_sum = 0;
        while (true) {
            size_t idx = task_index.fetch_add(1, memory_order_relaxed);
            if (idx >= tasks.size()) break;

            const Task& task = tasks[idx];
            const Pattern& pat = patterns[task.pattern_idx];
            const ValList& a_list = *pat.a_vals;
            const ValList& b_list = *pat.b_vals;
            const ValList& c_list = *pat.c_vals;
            const ValList& d_list = *pat.d_vals;

            int g = pat.g;
            int nmax = pat.nmax;

            for (int ai = task.a_start; ai < task.a_end; ++ai) {
                int a = a_list.v[ai];
                int a2 = a_list.v2[ai];
                bool a_zero = (a == 0);

                for (size_t bi = 0; bi < b_list.v.size(); ++bi) {
                    int b = b_list.v[bi];
                    if (a_zero && b < 0) continue;
                    int b2 = b_list.v2[bi];
                    int s1 = a2 + b2;
                    if (s1 > nmax) break;

                    int gcd_ab = std::gcd(iabs(a), iabs(b));
                    bool ab_zero = a_zero && (b == 0);

                    for (size_t ci = 0; ci < c_list.v.size(); ++ci) {
                        int c = c_list.v[ci];
                        if (ab_zero && c < 0) continue;
                        int c2 = c_list.v2[ci];
                        int s2 = s1 + c2;
                        if (s2 > nmax) break;

                        int gcd_abc = std::gcd(gcd_ab, (int)iabs(c));
                        bool abc_zero = ab_zero && (c == 0);

                        for (size_t di = 0; di < d_list.v.size(); ++di) {
                            int d = d_list.v[di];
                            if (abc_zero && d <= 0) continue;
                            int d2 = d_list.v2[di];
                            int N = s2 + d2;
                            if (N > nmax) break;
                            if (N == 0) continue;

                            if (std::gcd(gcd_abc, (int)iabs(d)) != 1) continue;

                            int64_t u1 = (int64_t)(a2 + b2 - c2 - d2) / g;
                            int64_t u2 = (int64_t)(2LL * (b * c - a * d)) / g;
                            int64_t u3 = (int64_t)(2LL * (b * d + a * c)) / g;

                            int64_t v1 = (int64_t)(2LL * (b * c + a * d)) / g;
                            int64_t v2 = (int64_t)(a2 - b2 + c2 - d2) / g;
                            int64_t v3 = (int64_t)(2LL * (c * d - a * b)) / g;

                            int64_t w1 = (int64_t)(2LL * (b * d - a * c)) / g;
                            int64_t w2 = (int64_t)(2LL * (c * d + a * b)) / g;
                            int64_t w3 = (int64_t)(a2 - b2 - c2 + d2) / g;

                            // Axis span equals the sum of absolute components for that axis.
                            int wx = (int)(iabs(u1) + iabs(v1) + iabs(w1));
                            int wy = (int)(iabs(u2) + iabs(v2) + iabs(w2));
                            int wz = (int)(iabs(u3) + iabs(v3) + iabs(w3));

                            if (wx > n || wy > n || wz > n) continue;
                            int maxw = max(wx, max(wy, wz));
                            int smax = n / maxw;
                            if (smax <= 0) continue;

                            int64_t L0 = (int64_t)N / g;
                            int gu = gcd3(u1, u2, u3);
                            int gv = gcd3(v1, v2, v3);
                            int gw = gcd3(w1, w2, w3);
                            int64_t G0 = (int64_t)gu + gv + gw;

                            int64_t ax = (int64_t)n + 1 - wx;
                            int64_t ay = (int64_t)n + 1 - wy;
                            int64_t az = (int64_t)n + 1 - wz;
                            int64_t L = L0;
                            int64_t G = G0;

                            for (int s = 1; s <= smax; ++s) {
                                unsigned long long placements = (unsigned long long)ax * (unsigned long long)ay *
                                                               (unsigned long long)az;

                                unsigned long long L64 = (unsigned long long)L;
                                unsigned long long G64 = (unsigned long long)G;
                                // Lattice points in a lattice cube: L^3 + (L+1) * (gu+gv+gw) + 1.
                                unsigned long long points = L64 * L64 * L64 + (L64 + 1ULL) * G64 + 1ULL;

                                local_sum += (unsigned __int128)placements * (unsigned __int128)points;

                                ax -= wx;
                                ay -= wy;
                                az -= wz;
                                L += L0;
                                G += G0;
                            }
                        }
                    }
                }
            }
        }
        sums[tid] = local_sum;
    };

    vector<thread> pool;
    pool.reserve(threads);
    for (unsigned int t = 0; t < threads; ++t) {
        pool.emplace_back(worker, t);
    }
    for (auto& th : pool) th.join();

    unsigned __int128 total = 0;
    for (unsigned int t = 0; t < threads; ++t) total += sums[t];
    return total;
}

unsigned __int128 compute_S_exact(int n, unsigned int threads) {
    unsigned __int128 sum_frames = compute_sum_frames(n, threads);
    if (sum_frames % 24 != 0) {
        cerr << "Internal error: sum not divisible by 24.\n";
        exit(1);
    }
    return sum_frames / 24;
}

}  // namespace

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    auto check = [&](int n, unsigned long long expected) {
        unsigned __int128 got = compute_S_exact(n, 1);
        unsigned long long got64 = (unsigned long long)got;
        if (got64 != expected) {
            cerr << "Validation failed: S(" << n << ") got " << got64
                 << " expected " << expected << "\n";
            exit(1);
        }
    };

    check(1, 8ULL);
    check(2, 91ULL);
    check(4, 1878ULL);
    check(5, 5832ULL);
    check(10, 387003ULL);
    check(50, 29948928129ULL);
    cerr << "Validation checkpoints passed.\n";

    const int n = 5000;
    unsigned int threads = thread::hardware_concurrency();
    if (threads == 0) threads = 4;

    unsigned __int128 total = compute_S_exact(n, threads);
    unsigned long long mod = (unsigned long long)(total % 1000000000ULL);
    cout << mod << "\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.nio.file.*;
import java.util.*;
import java.util.regex.*;

public class Euler579 {
    private static final Pattern ANSWER_RE = Pattern.compile("answer\\s*:\\s*(.+)$", Pattern.CASE_INSENSITIVE);
    private static final Pattern EQUAL_RE = Pattern.compile("=\\s*(.+)$");

    private static String parseOutput(String stdout) {
        String[] lines = stdout.split("\\R");
        List<String> nonEmpty = new ArrayList<>();
        for (String line : lines) {
            String t = line.trim();
            if (!t.isEmpty()) {
                nonEmpty.add(t);
            }
        }
        if (nonEmpty.isEmpty()) {
            return "";
        }

        List<String> answers = new ArrayList<>();
        List<String> equals = new ArrayList<>();
        for (String line : nonEmpty) {
            Matcher m1 = ANSWER_RE.matcher(line);
            if (m1.find()) {
                answers.add(m1.group(1).trim());
            }
            Matcher m2 = EQUAL_RE.matcher(line);
            if (m2.find()) {
                equals.add(m2.group(1).trim());
            }
        }

        if (!answers.isEmpty()) {
            return answers.get(answers.size() - 1);
        }
        if (!equals.isEmpty()) {
            return equals.get(equals.size() - 1);
        }
        return nonEmpty.get(nonEmpty.size() - 1);
    }

    private static String pickCompiler() throws Exception {
        for (String compiler : List.of("clang++", "g++")) {
            Process probe = new ProcessBuilder("bash", "-lc", "command -v " + compiler)
                    .redirectErrorStream(true)
                    .start();
            String out = new String(probe.getInputStream().readAllBytes());
            int rc = probe.waitFor();
            if (rc == 0 && !out.trim().isEmpty()) {
                return compiler;
            }
        }
        throw new RuntimeException("No C++ compiler found (clang++/g++).");
    }

    private static Path cppSource(Path root) {
        return root.resolve("solutionsCpp").resolve("Euler579.cpp");
    }

    private static boolean shouldSkipCheckpoints(Path root) {
        Path src = cppSource(root);
        try {
            String text = Files.readString(src);
            return text.contains("--skip-checkpoints");
        } catch (Exception ex) {
            return false;
        }
    }

    private static Path ensureBridgeBinary() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = root.resolve("solutionsCpp").resolve(".euler579_java_bridge");

        boolean rebuild = Files.notExists(bin)
                || Files.getLastModifiedTime(src).compareTo(Files.getLastModifiedTime(bin)) > 0;

        if (rebuild) {
            String compiler = pickCompiler();
            Process compile = new ProcessBuilder(
                    compiler,
                    "-std=c++17",
                    "-O2",
                    src.toString(),
                    "-o",
                    bin.toString())
                    .inheritIO()
                    .start();
            if (compile.waitFor() != 0) {
                throw new RuntimeException("Failed to compile Euler579 C++ bridge.");
            }
        }

        return bin;
    }

    private static String runBridge(Path bin, Path root, Path srcDir) throws Exception {
        List<String> cmd = new ArrayList<>();
        cmd.add(bin.toString());
        if (shouldSkipCheckpoints(root)) {
            cmd.add("--skip-checkpoints");
        }

        Process first = new ProcessBuilder(cmd)
                .directory(root.toFile())
                .redirectErrorStream(true)
                .start();
        String out = new String(first.getInputStream().readAllBytes());
        int rc = first.waitFor();
        if (rc == 0) {
            return out;
        }

        Process second = new ProcessBuilder(cmd)
                .directory(srcDir.toFile())
                .redirectErrorStream(true)
                .start();
        String out2 = new String(second.getInputStream().readAllBytes());
        int rc2 = second.waitFor();
        if (rc2 == 0) {
            return out2;
        }

        throw new RuntimeException("Euler579 C++ bridge failed.\n" + out + "\n" + out2);
    }

    private static String solveViaCppBridge() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = ensureBridgeBinary();
        String out = runBridge(bin, root, src.getParent());
        String parsed = parseOutput(out);
        if (parsed.isEmpty()) {
            throw new RuntimeException("Euler579 C++ bridge produced empty output.");
        }
        return parsed;
    }

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