Problem 578: Integers with Decreasing Prime Powers

View on Project Euler

Project Euler Problem 578 Solution

EulerSolve provides an optimized solution for Project Euler Problem 578, Integers with Decreasing Prime Powers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Write a positive integer \(n\gt 1\) in the form $$n=\prod_{i=1}^{k} p_i^{e_i},\qquad p_1 \lt p_2 \lt \cdots \lt p_k,$$ where the primes are listed in increasing order. Problem 578 asks for the counting function $$C(N)=\#\left\{n\le N : e_1\ge e_2\ge \cdots \ge e_k\ge 1\right\},$$ with \(1\) also counted. The difficulty is that \(N\) is enormous, so the solution cannot test each integer separately. Instead, it counts valid factorizations by their exponent pattern and then counts how many prime choices realize each pattern. Mathematical Approach The implementations split the problem into two layers: first enumerate every feasible nonincreasing exponent pattern, then count the prime tuples that fit that pattern under the bound \(N\). Step 1: Separate Numbers by Exponent Signature Every valid integer \(n\gt 1\) has a unique signature $$\lambda=(e_1,e_2,\dots,e_k),\qquad e_1\ge e_2\ge \cdots \ge e_k\ge 1.$$ For a fixed signature define $$A_\lambda(N)=\#\left\{(p_1,\dots,p_k): p_1\lt \cdots \lt p_k,\ \prod_{i=1}^{k} p_i^{e_i}\le N\right\}.$$ Then $$C(N)=1+\sum_{\lambda} A_\lambda(N),$$ where the \(+1\) accounts for \(n=1\). A signature is feasible only if even its smallest possible realization fits: $$2^{e_1}3^{e_2}5^{e_3}\cdots p_k^{e_k}\le N.$$ This observation is the foundation for the signature-generation recursion....

Detailed mathematical approach

Problem Summary

Write a positive integer \(n\gt 1\) in the form

$$n=\prod_{i=1}^{k} p_i^{e_i},\qquad p_1 \lt p_2 \lt \cdots \lt p_k,$$

where the primes are listed in increasing order. Problem 578 asks for the counting function

$$C(N)=\#\left\{n\le N : e_1\ge e_2\ge \cdots \ge e_k\ge 1\right\},$$

with \(1\) also counted. The difficulty is that \(N\) is enormous, so the solution cannot test each integer separately. Instead, it counts valid factorizations by their exponent pattern and then counts how many prime choices realize each pattern.

Mathematical Approach

The implementations split the problem into two layers: first enumerate every feasible nonincreasing exponent pattern, then count the prime tuples that fit that pattern under the bound \(N\).

Step 1: Separate Numbers by Exponent Signature

Every valid integer \(n\gt 1\) has a unique signature

$$\lambda=(e_1,e_2,\dots,e_k),\qquad e_1\ge e_2\ge \cdots \ge e_k\ge 1.$$

For a fixed signature define

$$A_\lambda(N)=\#\left\{(p_1,\dots,p_k): p_1\lt \cdots \lt p_k,\ \prod_{i=1}^{k} p_i^{e_i}\le N\right\}.$$

Then

$$C(N)=1+\sum_{\lambda} A_\lambda(N),$$

where the \(+1\) accounts for \(n=1\). A signature is feasible only if even its smallest possible realization fits:

$$2^{e_1}3^{e_2}5^{e_3}\cdots p_k^{e_k}\le N.$$

This observation is the foundation for the signature-generation recursion.

Step 2: Count Prime Choices for One Fixed Signature

Let \(p_1,p_2,\dots\) denote the prime numbers. Suppose we have already committed to primes up through \(p_i\), so the next chosen prime must be larger than \(p_i\). For a remaining exponent list \((f_1,\dots,f_r)\) and a remaining limit \(L\), define

$$R(f_1,\dots,f_r;L,i)$$

to be the number of strictly increasing prime choices \(q_1\lt \cdots \lt q_r\), all greater than \(p_i\), such that

$$q_1^{f_1}q_2^{f_2}\cdots q_r^{f_r}\le L.$$

The recursive relation is

$$R(f_1,\dots,f_r;L,i)=\sum_{j>i,\ p_j^{f_1}\le L} R(f_2,\dots,f_r;L/p_j^{f_1},j).$$

The desired count for a full signature is simply

$$A_{(e_1,\dots,e_k)}(N)=R(e_1,\dots,e_k;N,0).$$

Step 3: Collapse Terminal Branches with \(\pi(x)\)

When only one exponent \(e\) remains, the recursion reduces to counting primes:

$$R(e;L,i)=\max\left(0,\pi\left(\left\lfloor L^{1/e}\right\rfloor\right)-i\right).$$

The reason is simple: any prime \(q\) larger than the first \(i\) primes is valid exactly when \(q^e\le L\).

When two exponents \(e\ge f\) remain, we obtain

$$R(e,f;L,i)=\sum_{j>i,\ p_j^e\le L}\max\left(0,\pi\left(\left\lfloor (L/p_j^e)^{1/f}\right\rfloor\right)-j\right).$$

This is much faster than descending one level deeper for every pair. The C++ implementation also gives special treatment to common exponents \(1,2,4,8\), because the corresponding bounds can be obtained by exact integer roots.

Step 4: Prune by the Minimal Possible Continuation

For long signatures, most branches are impossible. If we are considering the next prime \(p_j\) for the remaining exponents \((f_1,\dots,f_r)\), then even the smallest possible continuation would use consecutive larger primes:

$$p_j^{f_1}p_{j+1}^{f_2}\cdots p_{j+r-1}^{f_r}.$$

If this minimal product already exceeds \(L\), then no later choice of \(p_j\) can work, because later primes are only larger. So the search can stop immediately at that point instead of exploring dead branches.

Step 5: Enumerate All Feasible Signatures Exactly Once

A separate recursion generates the signatures themselves. Start with an exponent \(e_1\ge 1\) satisfying \(2^{e_1}\le N\). For the next position choose \(e_2\) with

$$1\le e_2\le e_1,\qquad 3^{e_2}\le N/2^{e_1}.$$

Then continue with \(5,7,11,\dots\) as the smallest possible future primes. Every time a prefix \((e_1,\dots,e_t)\) is feasible with these minimal primes, that prefix represents a genuine signature, so its contribution \(A_{(e_1,\dots,e_t)}(N)\) is added immediately. Extending the prefix by a new exponent \(e_{t+1}\le e_t\) guarantees that each nonincreasing signature appears once and only once.

Worked Example: Signature \((2,1)\) for \(N=100\)

Here we count numbers of the form \(p^2q\) with \(p\lt q\) and \(p^2q\le 100\). The formula becomes

$$A_{(2,1)}(100)=\sum_{p^2\le 100}\#\left\{q>p:\ q\le 100/p^2,\ q\text{ prime}\right\}.$$

Now evaluate each possible first prime:

$$\begin{aligned} p=2&:\quad q\le 25 &&\Rightarrow 8,\\ p=3&:\quad q\le 11 &&\Rightarrow 3,\\ p\ge 5&:\quad p^2q>100 &&\Rightarrow 0. \end{aligned}$$

So \(A_{(2,1)}(100)=11\). The same pruning idea also shows immediately that the longer signature \((2,2,1)\) is impossible at \(N=100\), because its minimal realization would be

$$2^2\cdot 3^2\cdot 5=180>100.$$

When all feasible signatures are summed, the total is \(C(100)=94\), which matches the implementation's checkpoint.

How the Code Works

The C++ implementation begins by precomputing a compressed prime-count structure. It stores enough information to answer \(\pi(x)\) quickly both for small values of \(x\) and for quotient values of the form \(\lfloor N/v\rfloor\). It also precomputes prime powers \(p^e\le N\) for the primes that can matter in recursive branches.

After that, one recursive routine generates every feasible nonincreasing exponent signature, using the minimal primes \(2,3,5,\dots\) only as a feasibility test. Another recursive routine counts prime tuples for the current signature. Whenever only one or two exponents remain, the search switches to direct prime-count formulas instead of recursing blindly. Independent top-level branches, distinguished by the first exponent, are processed in parallel and then added together.

The Python and Java implementations do not reimplement the mathematics separately. They invoke the same compiled computation and extract the final numeric answer, so all three language versions use the same counting method.

Complexity Analysis

There is no clean closed form depending only on \(N\), because the search cost is governed by how many exponent signatures survive the feasibility tests and how many prime branches survive pruning. The prime-count preprocessing stores \(O(\sqrt{N})\) distinct arguments and uses \(O(\sqrt{N})\) memory; its arithmetic cost is close to \(O(\sqrt{N}\log\log N)\).

After preprocessing, let \(B(N)\) denote the total number of recursive states that remain after pruning. Then the overall running time is roughly

$$O\!\left(\sqrt{N}\log\log N + B(N)\right),$$

with terminal branches answered by fast \(\pi(x)\) queries rather than full subtree expansion. In practice, this is dramatically smaller than enumerating all integers up to \(N\), which is why the method is fast enough for the Project Euler bound.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=578
  2. Prime-counting function: Wikipedia - Prime-counting function
  3. Fundamental theorem of arithmetic: Wikipedia - Fundamental theorem of arithmetic
  4. Integer partition: Wikipedia - Partition (number theory)
  5. Prime power: Wikipedia - Prime power

Problem 578 source code

C++

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

namespace {

using u64 = std::uint64_t;

u64 isqrt_u64(u64 x) {
    u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(x)));
    while ((r + 1) <= x / (r + 1)) ++r;
    while (r > x / r) --r;
    return r;
}

class Solver578 {
public:
    explicit Solver578(u64 n) : n_(n), r_(isqrt_u64(n)), b_(n / r_) {
        std::vector<u64> values;
        values.reserve(static_cast<std::size_t>(r_ + b_));
        for (u64 v = 1; v <= r_; ++v) values.push_back(n_ / v);
        for (u64 v = b_; v >= 1; --v) {
            if (v == b_) continue;
            values.push_back(v);
        }

        count0_.resize(static_cast<std::size_t>(b_));
        for (u64 v = 1; v <= b_; ++v) count0_[static_cast<std::size_t>(v - 1)] = v - 1;

        count1_.resize(static_cast<std::size_t>(r_));
        for (u64 v = 1; v <= r_; ++v) count1_[static_cast<std::size_t>(v - 1)] = n_ / v - 1;

        for (u64 p = 2; p <= r_; ++p) {
            const u64 prev = count0_[static_cast<std::size_t>(p - 2)];
            if (count0_[static_cast<std::size_t>(p - 1)] <= prev) continue;

            primes_.push_back(p);
            const u64 p2 = p * p;

            for (u64 v : values) {
                if (v < p2) break;

                if (v <= b_) {
                    count0_[static_cast<std::size_t>(v - 1)] -=
                        count0_[static_cast<std::size_t>(v / p - 1)] - prev;
                } else {
                    const std::size_t idx = static_cast<std::size_t>(n_ / v - 1);
                    const u64 vp = v / p;
                    if (vp <= b_) {
                        count1_[idx] -= count0_[static_cast<std::size_t>(vp - 1)] - prev;
                    } else {
                        count1_[idx] -= count1_[static_cast<std::size_t>(n_ / vp - 1)] - prev;
                    }
                }
            }
        }

        prime_pows_.resize(primes_.size());
        for (std::size_t idx = 0; idx < primes_.size(); ++idx) {
            const u64 p = primes_[idx];
            auto& pw = prime_pows_[idx];
            pw.reserve(8);
            pw.push_back(1);
            while (true) {
                const u64 cur = pw.back();
                if (cur > n_ / p) break;
                pw.push_back(cur * p);
            }
        }
    }

    u64 solve(unsigned int thread_count = 0) const {
        if (primes_.empty()) return 1;
        std::vector<std::pair<unsigned char, u64>> roots;
        roots.reserve(64);
        const u64 first_prime = primes_[0];
        u64 p_pow = first_prime;
        for (int e = 1; e <= 100 && p_pow <= n_; ++e) {
            roots.emplace_back(static_cast<unsigned char>(e), p_pow);
            if (p_pow > n_ / first_prime) break;
            p_pow *= first_prime;
        }

        if (roots.empty()) return 1;

        if (thread_count == 0) {
            thread_count = std::thread::hardware_concurrency();
            if (thread_count == 0) thread_count = 1;
        }
        thread_count = std::min<unsigned int>(thread_count, static_cast<unsigned int>(roots.size()));

        if (thread_count <= 1) {
            return 1 + solve_roots_sequential(roots);
        }

        std::atomic<std::size_t> next(0);
        std::vector<u64> partial(static_cast<std::size_t>(thread_count), 0);
        std::vector<std::thread> workers;
        workers.reserve(static_cast<std::size_t>(thread_count));

        for (unsigned int t = 0; t < thread_count; ++t) {
            workers.emplace_back([&, t]() {
                std::vector<unsigned char> current;
                current.reserve(64);
                u64 local = 0;
                while (true) {
                    const std::size_t idx = next.fetch_add(1, std::memory_order_relaxed);
                    if (idx >= roots.size()) break;
                    const auto [e, pw] = roots[idx];
                    current.clear();
                    current.push_back(e);
                    local += signature_count(current, 0, n_, 0);
                    local += enumerate_sum(n_ / pw, 1, e, current);
                }
                partial[static_cast<std::size_t>(t)] = local;
            });
        }
        for (auto& th : workers) th.join();
        u64 total = 1;
        for (u64 v : partial) total += v;
        return total;
    }

private:
    u64 n_;
    u64 r_;
    u64 b_;
    std::vector<u64> primes_;
    std::vector<u64> count0_;
    std::vector<u64> count1_;
    std::vector<std::vector<u64>> prime_pows_;

    bool pow_leq_idx(std::size_t idx, int e, u64 limit, u64& out) const {
        if (e < 0) return false;
        const auto& pw = prime_pows_[idx];
        if (static_cast<std::size_t>(e) >= pw.size()) return false;
        const u64 val = pw[static_cast<std::size_t>(e)];
        if (val > limit) return false;
        out = val;
        return true;
    }

    u64 pi(u64 v) const {
        if (v <= 1) return 0;
        if (v <= b_) return count0_[static_cast<std::size_t>(v - 1)];
        return count1_[static_cast<std::size_t>(n_ / v - 1)];
    }

    u64 count_pairs(int e0, int e1, u64 limit, std::size_t i) const {
        if (i + 1 >= primes_.size()) return 0;

        u64 total = 0;
        for (std::size_t j = i; j < primes_.size(); ++j) {
            u64 p0e = 1;
            if (!pow_leq_idx(j, e0, limit, p0e)) break;

            const u64 lim1 = limit / p0e;
            if (e1 == 1) {
                const u64 cnt = pi(lim1);
                if (cnt <= j + 1) break;
                total += cnt - (j + 1);
            } else if (e1 == 2) {
                const u64 cnt = pi(isqrt_u64(lim1));
                if (cnt <= j + 1) break;
                total += cnt - (j + 1);
            } else if (e1 == 4) {
                const u64 cnt = pi(isqrt_u64(isqrt_u64(lim1)));
                if (cnt <= j + 1) break;
                total += cnt - (j + 1);
            } else if (e1 == 8) {
                const u64 cnt = pi(isqrt_u64(isqrt_u64(isqrt_u64(lim1))));
                if (cnt <= j + 1) break;
                total += cnt - (j + 1);
            } else {
                for (std::size_t k = j + 1; k < primes_.size(); ++k) {
                    u64 p1e = 1;
                    if (!pow_leq_idx(k, e1, lim1, p1e)) break;
                    ++total;
                }
            }
        }

        return total;
    }

    u64 signature_count(const std::vector<unsigned char>& s, int pos, u64 limit,
                        std::size_t i) const {
        const int rem = static_cast<int>(s.size()) - pos;
        if (rem == 1) {
            const int e = static_cast<int>(s[pos]);
            if (e == 1) {
                const u64 cnt = pi(limit);
                return (cnt > i) ? cnt - i : 0ULL;
            }
            if (e == 2) {
                const u64 cnt = pi(isqrt_u64(limit));
                return (cnt > i) ? cnt - i : 0ULL;
            }
            if (e == 4) {
                const u64 cnt = pi(isqrt_u64(isqrt_u64(limit)));
                return (cnt > i) ? cnt - i : 0ULL;
            }
            if (e == 8) {
                const u64 cnt = pi(isqrt_u64(isqrt_u64(isqrt_u64(limit))));
                return (cnt > i) ? cnt - i : 0ULL;
            }

            u64 cnt = 0;
            for (std::size_t j = i; j < primes_.size(); ++j) {
                u64 pw = 1;
                if (!pow_leq_idx(j, e, limit, pw)) break;
                ++cnt;
            }
            return cnt;
        }
        if (rem == 2) {
            return count_pairs(static_cast<int>(s[pos]), static_cast<int>(s[pos + 1]), limit, i);
        }

        u64 total = 0;
        for (std::size_t j = i; j < primes_.size(); ++j) {
            bool ok = true;
            u64 minimal = 1;
            const std::size_t max_t =
                std::min<std::size_t>(static_cast<std::size_t>(rem), primes_.size() - j);

            for (std::size_t t = 0; t < max_t; ++t) {
                u64 pw = 1;
                if (!pow_leq_idx(j + t, static_cast<int>(s[pos + static_cast<int>(t)]),
                                 limit / minimal, pw)) {
                    ok = false;
                    break;
                }
                minimal *= pw;
            }
            if (!ok) break;

            u64 first_pw = 1;
            if (!pow_leq_idx(j, static_cast<int>(s[pos]), limit, first_pw)) break;

            total += signature_count(s, pos + 1, limit / first_pw, j + 1);
        }

        return total;
    }

    u64 solve_roots_sequential(const std::vector<std::pair<unsigned char, u64>>& roots) const {
        std::vector<unsigned char> current;
        current.reserve(64);
        u64 total = 0;
        for (const auto [e, pw] : roots) {
            current.clear();
            current.push_back(e);
            total += signature_count(current, 0, n_, 0);
            total += enumerate_sum(n_ / pw, 1, e, current);
        }
        return total;
    }

    u64 enumerate_sum(u64 limit, std::size_t i, int exponent_bound,
                      std::vector<unsigned char>& current) const {
        if (i >= primes_.size()) return 0;
        u64 total = 0;

        const u64 p = primes_[i];
        u64 p_pow = p;

        for (int e = 1; e <= exponent_bound && p_pow <= limit; ++e) {
            current.push_back(static_cast<unsigned char>(e));
            total += signature_count(current, 0, n_, 0);
            total += enumerate_sum(limit / p_pow, i + 1, e, current);
            current.pop_back();

            if (p_pow > limit / p) break;
            p_pow *= p;
        }
        return total;
    }
};

bool is_decreasing_prime_power(u64 n) {
    if (n <= 1) return true;
    std::vector<int> exponents;
    u64 x = n;

    for (u64 p = 2; p * p <= x; ++p) {
        if (x % p != 0) continue;
        int cnt = 0;
        while (x % p == 0) {
            x /= p;
            ++cnt;
        }
        exponents.push_back(cnt);
    }

    if (x > 1) exponents.push_back(1);

    for (std::size_t i = 1; i < exponents.size(); ++i) {
        if (exponents[i - 1] < exponents[i]) return false;
    }
    return true;
}

u64 brute_count(u64 n) {
    u64 cnt = 0;
    for (u64 i = 1; i <= n; ++i) {
        if (is_decreasing_prime_power(i)) ++cnt;
    }
    return cnt;
}

bool run_validations() {
    struct Check {
        u64 n;
        u64 expected;
    };

    const std::vector<Check> checks = {
        {100ULL, 94ULL},
        {1'000'000ULL, 922'052ULL},
    };

    for (const auto& chk : checks) {
        Solver578 solver(chk.n);
        const u64 got = solver.solve();
        if (got != chk.expected) {
            std::cerr << "Validation failed for C(" << chk.n << "): got " << got
                      << ", expected " << chk.expected << "\n";
            return false;
        }
    }

    const u64 small_n = 10'000ULL;
    Solver578 small_solver(small_n);
    const u64 got_small = small_solver.solve();
    const u64 expected_small = brute_count(small_n);
    if (got_small != expected_small) {
        std::cerr << "Validation failed for brute-force C(" << small_n << "): got "
                  << got_small << ", expected " << expected_small << "\n";
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    u64 n = 10'000'000'000'000ULL;
    if (argc > 1) n = std::strtoull(argv[1], nullptr, 10);

    if (!run_validations()) return 1;

    Solver578 solver(n);
    std::cout << solver.solve() << '\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 Euler578 {
    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("Euler578.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(".euler578_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 Euler578 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("Euler578 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("Euler578 C++ bridge produced empty output.");
        }
        return parsed;
    }

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