Problem 667: Moving Pentagon

Project Euler Problem 667 Solution

EulerSolve provides an optimized solution for Project Euler Problem 667, Moving Pentagon, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We study equilateral pentagons with side length $1$ that must move through a right-angled corridor of width $W$. For such a pentagon let $A$ be its area, and let $W$ be the smallest corridor width that still allows a continuous motion from one straight branch of the corridor to the other. The objective is to maximize $$\frac{A}{W^2}.$$ This ratio is scale invariant: scaling the pentagon by a factor $\lambda$ multiplies the area by $\lambda^2$ and the required corridor width by $\lambda$. So it is enough to work with unit-edge pentagons and optimize a normalized shape parameter. Mathematical Approach The implementation solves the problem numerically but the geometry is highly structured. First it restricts to a symmetric one-parameter family of equilateral pentagons, then it computes the minimal corridor width for each candidate by sampling all orientations and testing whether those orientations form a continuous feasible passage around the corner....

Detailed mathematical approach

Problem Summary

We study equilateral pentagons with side length $1$ that must move through a right-angled corridor of width $W$. For such a pentagon let $A$ be its area, and let $W$ be the smallest corridor width that still allows a continuous motion from one straight branch of the corridor to the other.

The objective is to maximize

$$\frac{A}{W^2}.$$

This ratio is scale invariant: scaling the pentagon by a factor $\lambda$ multiplies the area by $\lambda^2$ and the required corridor width by $\lambda$. So it is enough to work with unit-edge pentagons and optimize a normalized shape parameter.

Mathematical Approach

The implementation solves the problem numerically but the geometry is highly structured. First it restricts to a symmetric one-parameter family of equilateral pentagons, then it computes the minimal corridor width for each candidate by sampling all orientations and testing whether those orientations form a continuous feasible passage around the corner.

Step 1: Reduce the pentagon to one angle

Symmetry about the perpendicular bisector of the first edge lets us place the first four visible vertices as

$$v_0=(0,0),\qquad v_1=(1,0),\qquad v_2=(1+\cos a,\sin a),\qquad v_4=(-\cos a,\sin a).$$

This already guarantees

$$|v_1-v_0|=|v_2-v_1|=|v_0-v_4|=1.$$

The only missing point is $v_3$, which must satisfy

$$|v_3-v_2|=|v_3-v_4|=1.$$

So the entire symmetric equilateral pentagon is controlled by the single angle parameter $a$.

Step 2: Recover the fifth vertex from two unit circles

The distance between the two circle centers is

$$d=|v_2-v_4|=1+2\cos a.$$

For the two unit circles to intersect we need $d\le 2$, hence this construction requires

$$a\ge \frac{\pi}{3}.$$

The midpoint of $v_2v_4$ is $\left(\frac12,\sin a\right)$, so the two possible intersection points are

$$v_3=\left(\frac12,\sin a \pm \sqrt{1-\frac{d^2}{4}}\right).$$

These are the two mirror candidates examined by the implementation. For each candidate, the area is computed with the shoelace formula

$$A=\frac12\left|\sum_{i=0}^{4}\left(x_i y_{i+1}-x_{i+1} y_i\right)\right|,$$

with indices taken cyclically.

Step 3: Translate the corridor constraints into width profiles

Fix an orientation angle $\theta$, rotate the pentagon by $\theta$, and then translate the rotated coordinates so that the minimum $x$-coordinate and minimum $y$-coordinate are both $0$. For the translated polygon $P_\theta$, define

$$w_x(\theta)=x_{\max}-x_{\min},\qquad w_y(\theta)=y_{\max}-y_{\min}.$$

These are the widths needed for the pentagon to fit inside a straight vertical corridor and a straight horizontal corridor of width $W$, respectively.

The corner constraint is different. In the translated picture, a right-angled corridor of width $W$ occupies the union of the strips $x\le W$ and $y\le W$, so the forbidden region is the upper-right square where both coordinates exceed $W$. Therefore the pentagon fits at the corner exactly when

$$\max_{p\in P_\theta}\min(p_x,p_y)\le W.$$

This motivates the corner-width functional

$$w(\theta)=\max_{p\in P_\theta}\min(p_x,p_y).$$

Step 4: Why only vertices and diagonal crossings matter

Along any edge segment of the polygon, the function $\min(x,y)$ is piecewise linear. Away from the diagonal $x=y$, it is simply one coordinate or the other. The only place where the active branch can switch is where the edge crosses the diagonal.

So the maximum value of $\min(x,y)$ on the polygon boundary is attained either at a vertex or at a point where an edge meets the diagonal $x=y$.

The implementation uses exactly that fact: for every sampled orientation it checks all shifted vertices and all edge-diagonal intersections, and the largest such value becomes $w(\theta)$.

Step 5: Turn continuous motion into connectivity on the angle circle

For a fixed corridor width $W$, let

$$S_W=\{\theta\in[0,2\pi): w(\theta)\le W\}.$$

These are the orientations that can be placed at the corner. A continuous passage around the bend requires more than isolated feasible orientations. We need one connected arc of $S_W$ that contains

an orientation with $w_y(\theta)\le W$, so the pentagon can lie inside a straight horizontal branch, and

an orientation with $w_x(\theta)\le W$, so it can lie inside a straight vertical branch.

Because orientation is periodic, the angle domain is a circle. The implementation samples that circle densely, sorts the samples by increasing $w(\theta)$, activates them in threshold order, and merges neighboring active samples. Each connected component stores the smallest observed $w_x$ and $w_y$. The first threshold whose component satisfies both minima is the required motion width for that pentagon.

Worked Example: A concrete symmetric pentagon at $a=\frac{\pi}{2}$

Take

$$a=\frac{\pi}{2}.$$

Then

$$v_0=(0,0),\qquad v_1=(1,0),\qquad v_2=(1,1),\qquad v_4=(0,1),$$

and the two circle intersections are

$$v_3=\left(\frac12,1\pm\frac{\sqrt3}{2}\right).$$

Using the upper intersection gives a simple pentagon with area

$$A=1+\frac{\sqrt3}{4}\approx 1.4330127.$$

At orientation $\theta=0$, the polygon is already translated with $x_{\min}=y_{\min}=0$, so

$$w_x(0)=1,\qquad w_y(0)=1+\frac{\sqrt3}{2},\qquad w(0)=1.$$

This means the pentagon can sit at the corner of a width-$1$ corridor, because no point has both coordinates larger than $1$. But it does not fit completely inside a straight horizontal branch of width $1$, since its vertical span exceeds $1$. This is exactly why the algorithm tracks all three profiles $w$, $w_x$, and $w_y$.

Step 6: Optimize the normalized objective

For each admissible $a$, the implementation evaluates both mirror candidates, computes the corresponding minimal corridor width $W(a)$, and keeps the larger value of

$$F(a)=\frac{A(a)}{W(a)^2}.$$

The search is one-dimensional, and within the narrow region containing the best symmetric pentagon the numerical profile is handled efficiently by a golden-section search.

How the Code Works

The C++ implementation precomputes a dense table of trigonometric values for equally spaced orientations around the full circle. For each candidate pentagon it rotates the five vertices at every sampled angle, shifts the pose so the lower-left touching position is normalized, and records the three geometric profiles $w$, $w_x$, and $w_y$.

It then sorts the sampled orientations by the corner-width profile $w$. As the threshold rises, neighboring active samples on the cyclic angle grid are merged with a disjoint-set structure. Each connected component remembers the smallest straight-corridor widths seen so far, so the first component with both $w_x\le W$ and $w_y\le W$ determines the motion width.

At the outer level, the implementation performs a golden-section search over the single angle parameter and evaluates both circle-intersection branches at each step. The final candidate is checked numerically to confirm that all five edges still have unit length and that the pentagon is simple. The Python and Java implementations delegate to the same compiled numerical core and return the parsed numeric result.

Complexity Analysis

Let $N_\theta$ be the number of sampled orientations. For one pentagon, computing all rotated poses and the three width profiles costs $O(N_\theta)$ geometric work because the polygon has only five vertices and five edges. Sorting the orientations by $w$ costs

$$O(N_\theta\log N_\theta),$$

and the connectivity sweep with a disjoint-set structure is effectively linear, $O(N_\theta \alpha(N_\theta))$.

So one motion-width evaluation is dominated by $O(N_\theta\log N_\theta)$ time and $O(N_\theta)$ memory. If the outer search uses $T$ objective evaluations, the total cost is

$$O(T\,N_\theta\log N_\theta).$$

The C++ version parallelizes the per-angle geometry across hardware threads, which improves wall-clock time but does not change the asymptotic bound.

Footnotes and References

Problem page: https://projecteuler.net/problem=667
Shoelace formula: Wikipedia - Shoelace formula
Moving sofa problem: Wikipedia - Moving sofa problem
Golden-section search: Wikipedia - Golden-section search
Disjoint-set data structure: Wikipedia - Disjoint-set data structure

Problem 667 source code

C++

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <thread>
#include <vector>

namespace {

struct Point {
    long double x;
    long double y;
};

constexpr long double kPi = 3.141592653589793238462643383279502884L;
constexpr long double kEps = 1e-13L;

long double dist(const Point& a, const Point& b) {
    long double dx = a.x - b.x;
    long double dy = a.y - b.y;
    return std::sqrt(dx * dx + dy * dy);
}

long double polygon_area(const std::vector<Point>& poly) {
    long double s = 0.0L;
    for (std::size_t i = 0; i < poly.size(); ++i) {
        const Point& p = poly[i];
        const Point& q = poly[(i + 1) % poly.size()];
        s += p.x * q.y - q.x * p.y;
    }
    return std::fabsl(s) * 0.5L;
}

long double orient(const Point& a, const Point& b, const Point& c) {
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}

bool on_segment(const Point& a, const Point& b, const Point& c) {
    return std::min(a.x, b.x) - kEps <= c.x && c.x <= std::max(a.x, b.x) + kEps &&
           std::min(a.y, b.y) - kEps <= c.y && c.y <= std::max(a.y, b.y) + kEps;
}

bool segments_intersect(const Point& a, const Point& b, const Point& c, const Point& d) {
    long double o1 = orient(a, b, c);
    long double o2 = orient(a, b, d);
    long double o3 = orient(c, d, a);
    long double o4 = orient(c, d, b);
    if (std::fabsl(o1) < kEps && on_segment(a, b, c)) return true;
    if (std::fabsl(o2) < kEps && on_segment(a, b, d)) return true;
    if (std::fabsl(o3) < kEps && on_segment(c, d, a)) return true;
    if (std::fabsl(o4) < kEps && on_segment(c, d, b)) return true;
    return (o1 * o2 < 0) && (o3 * o4 < 0);
}

bool is_simple_polygon(const std::vector<Point>& poly) {
    std::size_t n = poly.size();
    for (std::size_t i = 0; i < n; ++i) {
        Point a1 = poly[i];
        Point a2 = poly[(i + 1) % n];
        for (std::size_t j = i + 1; j < n; ++j) {
            if (j == i || (j + 1) % n == i || j == (i + 1) % n) {
                continue;
            }
            Point b1 = poly[j];
            Point b2 = poly[(j + 1) % n];
            if (segments_intersect(a1, a2, b1, b2)) {
                return false;
            }
        }
    }
    return true;
}

std::vector<std::vector<Point>> build_candidates(long double a) {
    long double b = kPi - a;
    Point v0{0.0L, 0.0L};
    Point v1{1.0L, 0.0L};
    Point v2{1.0L + std::cosl(a), std::sinl(a)};
    Point v4{std::cosl(b), std::sinl(b)};

    long double dx = v4.x - v2.x;
    long double dy = v4.y - v2.y;
    long double d = std::sqrt(dx * dx + dy * dy);
    long double half = 0.5L * d;
    long double h = std::sqrt(std::max(0.0L, 1.0L - half * half));
    long double mx = (v2.x + v4.x) * 0.5L;
    long double my = (v2.y + v4.y) * 0.5L;
    long double ox = -dy / d * h;
    long double oy = dx / d * h;

    Point p1{mx + ox, my + oy};
    Point p2{mx - ox, my - oy};
    std::vector<Point> poly1{v0, v1, v2, p1, v4};
    std::vector<Point> poly2{v0, v1, v2, p2, v4};
    return {poly1, poly2};
}

struct WidthData {
    std::vector<long double> w;
    std::vector<long double> wx;
    std::vector<long double> wy;
};

WidthData compute_widths(const std::vector<Point>& poly,
                         const std::vector<long double>& cos_t,
                         const std::vector<long double>& sin_t) {
    std::size_t n = cos_t.size();
    WidthData data;
    data.w.assign(n, 0.0L);
    data.wx.assign(n, 0.0L);
    data.wy.assign(n, 0.0L);

    std::size_t threads = std::max<unsigned>(1U, std::thread::hardware_concurrency());
    std::size_t chunk = (n + threads - 1) / threads;
    std::vector<std::thread> workers;
    workers.reserve(threads);

    auto worker = [&](std::size_t start) {
        std::size_t end = std::min(n, start + chunk);
        std::vector<Point> pts(poly.size());
        std::vector<Point> shifted(poly.size());
        for (std::size_t i = start; i < end; ++i) {
            long double c = cos_t[i];
            long double s = sin_t[i];
            long double minx = 1e100L, miny = 1e100L;
            long double maxx = -1e100L, maxy = -1e100L;
            for (std::size_t k = 0; k < poly.size(); ++k) {
                long double x = c * poly[k].x - s * poly[k].y;
                long double y = s * poly[k].x + c * poly[k].y;
                pts[k] = {x, y};
                minx = std::min(minx, x);
                miny = std::min(miny, y);
                maxx = std::max(maxx, x);
                maxy = std::max(maxy, y);
            }
            for (std::size_t k = 0; k < poly.size(); ++k) {
                shifted[k] = {pts[k].x - minx, pts[k].y - miny};
            }

            long double maxv = 0.0L;
            for (const auto& p : shifted) {
                long double v = (p.x < p.y) ? p.x : p.y;
                if (v > maxv) {
                    maxv = v;
                }
            }
            for (std::size_t k = 0; k < shifted.size(); ++k) {
                const Point& a = shifted[k];
                const Point& b = shifted[(k + 1) % shifted.size()];
                long double d1 = a.x - a.y;
                long double d2 = b.x - b.y;
                if (std::fabsl(d1) < kEps && std::fabsl(d2) < kEps) {
                    maxv = std::max(maxv, a.x);
                } else if (std::fabsl(d1) < kEps) {
                    maxv = std::max(maxv, a.x);
                } else if (std::fabsl(d2) < kEps) {
                    maxv = std::max(maxv, b.x);
                } else if (d1 * d2 < 0) {
                    long double t = d1 / (d1 - d2);
                    long double xi = a.x + t * (b.x - a.x);
                    maxv = std::max(maxv, xi);
                }
            }
            data.w[i] = maxv;
            data.wx[i] = maxx - minx;
            data.wy[i] = maxy - miny;
        }
    };

    for (std::size_t t = 0; t < threads; ++t) {
        std::size_t start = t * chunk;
        if (start >= n) break;
        workers.emplace_back(worker, start);
    }
    for (auto& th : workers) {
        th.join();
    }
    return data;
}

long double compute_min_width(const std::vector<Point>& poly,
                              const std::vector<long double>& cos_t,
                              const std::vector<long double>& sin_t) {
    WidthData data = compute_widths(poly, cos_t, sin_t);
    std::size_t n = data.w.size();

    std::vector<std::size_t> idx(n);
    std::iota(idx.begin(), idx.end(), 0);
    std::sort(idx.begin(), idx.end(),
              [&](std::size_t a, std::size_t b) { return data.w[a] < data.w[b]; });

    std::vector<std::size_t> parent(n);
    std::vector<bool> active(n, false);
    std::vector<long double> min_wx(n, 0.0L);
    std::vector<long double> min_wy(n, 0.0L);

    auto find = [&](auto self, std::size_t x) -> std::size_t {
        if (parent[x] == x) return x;
        parent[x] = self(self, parent[x]);
        return parent[x];
    };

    auto unite = [&](std::size_t a, std::size_t b) {
        std::size_t ra = find(find, a);
        std::size_t rb = find(find, b);
        if (ra == rb) return ra;
        parent[rb] = ra;
        min_wx[ra] = std::min(min_wx[ra], min_wx[rb]);
        min_wy[ra] = std::min(min_wy[ra], min_wy[rb]);
        return ra;
    };

    for (std::size_t id : idx) {
        active[id] = true;
        parent[id] = id;
        min_wx[id] = data.wx[id];
        min_wy[id] = data.wy[id];

        std::size_t left = (id + n - 1) % n;
        std::size_t right = (id + 1) % n;
        if (active[left]) {
            unite(id, left);
        }
        if (active[right]) {
            unite(id, right);
        }

        std::size_t root = find(find, id);
        long double w = data.w[id];
        if (min_wx[root] <= w + 1e-15L && min_wy[root] <= w + 1e-15L) {
            return w;
        }
    }
    return data.w[idx.back()];
}

struct EvalResult {
    long double value;
    long double area;
    long double width;
    std::vector<Point> poly;
};

EvalResult evaluate(long double a,
                    const std::vector<long double>& cos_t,
                    const std::vector<long double>& sin_t) {
    EvalResult best{0.0L, 0.0L, 0.0L, {}};
    for (const auto& poly : build_candidates(a)) {
        long double A = polygon_area(poly);
        long double W = compute_min_width(poly, cos_t, sin_t);
        long double val = A / (W * W);
        if (val > best.value) {
            best = {val, A, W, poly};
        }
    }
    return best;
}

long double objective(long double a,
                      const std::vector<long double>& cos_t,
                      const std::vector<long double>& sin_t) {
    return evaluate(a, cos_t, sin_t).value;
}

}  // namespace

int main() {
    const std::size_t N = 500000;
    std::vector<long double> cos_t(N), sin_t(N);
    for (std::size_t i = 0; i < N; ++i) {
        long double theta = 2.0L * kPi * static_cast<long double>(i) / static_cast<long double>(N);
        cos_t[i] = std::cosl(theta);
        sin_t[i] = std::sinl(theta);
    }

    // Symmetry around the mid-edge axis reduces the shape to one angle parameter.
    // Golden-section search for the optimal symmetric pentagon.
    long double lo = 1.055L;
    long double hi = 1.062L;
    const long double phi = (std::sqrt(5.0L) - 1.0L) * 0.5L;

    long double x1 = hi - phi * (hi - lo);
    long double x2 = lo + phi * (hi - lo);
    long double f1 = objective(x1, cos_t, sin_t);
    long double f2 = objective(x2, cos_t, sin_t);

    for (int iter = 0; iter < 60; ++iter) {
        if (f1 < f2) {
            lo = x1;
            x1 = x2;
            f1 = f2;
            x2 = lo + phi * (hi - lo);
            f2 = objective(x2, cos_t, sin_t);
        } else {
            hi = x2;
            x2 = x1;
            f2 = f1;
            x1 = hi - phi * (hi - lo);
            f1 = objective(x1, cos_t, sin_t);
        }
    }

    long double best_a = (lo + hi) * 0.5L;
    EvalResult result = evaluate(best_a, cos_t, sin_t);
    std::vector<Point> poly = result.poly;
    long double A = result.area;
    long double W = result.width;
    long double answer = result.value;

    // Validation checkpoints.
    long double min_len = 1e100L, max_len = 0.0L;
    for (std::size_t i = 0; i < poly.size(); ++i) {
        long double len = dist(poly[i], poly[(i + 1) % poly.size()]);
        min_len = std::min(min_len, len);
        max_len = std::max(max_len, len);
    }
    if (std::fabsl(max_len - min_len) > 5e-10L) {
        std::cerr << "Edge length mismatch: " << std::setprecision(18)
                  << min_len << " vs " << max_len << "\n";
        return 1;
    }
    if (!is_simple_polygon(poly)) {
        std::cerr << "Pentagon is not simple.\n";
        return 1;
    }
    std::cout.setf(std::ios::fixed);
    std::cout << std::setprecision(10) << static_cast<double>(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.nio.file.*;
import java.util.*;
import java.util.regex.*;

public class Euler667 {
    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("Euler667.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(".euler667_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 Euler667 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("Euler667 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("Euler667 C++ bridge produced empty output.");
        }
        return parsed;
    }

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