Problem 254: Sums of Digit Factorials

View on Project Euler

Project Euler Problem 254 Solution

EulerSolve provides an optimized solution for Project Euler Problem 254, Sums of Digit Factorials, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For a positive integer \(n\), define $$f(n)=\sum_{d\in\mathrm{digits}(n)} d!,$$ and then define $$sf(n)=s\bigl(f(n)\bigr),$$ where \(s(x)\) is the ordinary decimal digit sum. For each \(i\ge 1\), the problem introduces $$g(i)=\min\{n\ge 1: sf(n)=i\},\qquad sg(i)=s\bigl(g(i)\bigr).$$ The goal is to compute $$\sum_{i=1}^{150} sg(i).$$ A direct search over \(n\) is completely impractical. The code therefore changes viewpoint: it searches over the value $$y=f(n),$$ derives the smallest decimal \(n\) compatible with that \(y\), and uses a bitset dynamic program to find all relevant \(y\) grouped by decimal digit sum and residue modulo \(9!\). Mathematical Approach 1. Search Over \(y=f(n)\) Instead of Over \(n\) If the digit \(d\) appears \(c_d\) times in \(n\), then $$f(n)=\sum_{d=0}^{9} c_d\,d!.$$ So a number \(n\) is completely described, for the purposes of \(f\), by the multiplicities of its digits. Once a target value \(y\) is fixed, the question becomes: Among all digit multisets whose factorial sum is \(y\), which multiset yields the smallest decimal integer? For any fixed multiset of digits, the smallest corresponding integer is obtained by arranging the digits in nondecreasing order. So for each \(y\), we only need the lexicographically smallest valid multiset. 2....

Detailed mathematical approach

Problem Summary

For a positive integer \(n\), define

$$f(n)=\sum_{d\in\mathrm{digits}(n)} d!,$$

and then define

$$sf(n)=s\bigl(f(n)\bigr),$$

where \(s(x)\) is the ordinary decimal digit sum.

For each \(i\ge 1\), the problem introduces

$$g(i)=\min\{n\ge 1: sf(n)=i\},\qquad sg(i)=s\bigl(g(i)\bigr).$$

The goal is to compute

$$\sum_{i=1}^{150} sg(i).$$

A direct search over \(n\) is completely impractical. The code therefore changes viewpoint: it searches over the value

$$y=f(n),$$

derives the smallest decimal \(n\) compatible with that \(y\), and uses a bitset dynamic program to find all relevant \(y\) grouped by decimal digit sum and residue modulo \(9!\).

Mathematical Approach

1. Search Over \(y=f(n)\) Instead of Over \(n\)

If the digit \(d\) appears \(c_d\) times in \(n\), then

$$f(n)=\sum_{d=0}^{9} c_d\,d!.$$

So a number \(n\) is completely described, for the purposes of \(f\), by the multiplicities of its digits. Once a target value \(y\) is fixed, the question becomes:

Among all digit multisets whose factorial sum is \(y\), which multiset yields the smallest decimal integer?

For any fixed multiset of digits, the smallest corresponding integer is obtained by arranging the digits in nondecreasing order. So for each \(y\), we only need the lexicographically smallest valid multiset.

2. Why a Minimal Candidate Never Needs the Digit 0

This is one of the key structural facts behind the code.

First, note that

$$0!=1!=1.$$

If a candidate multiset contains both a 0 and a 1, then we may replace that pair by a single digit 2, because

$$0!+1!=1+1=2=2!.$$

This preserves \(f(n)\) while shortening the decimal length, so the original number cannot have been minimal.

If a candidate contains a 0 but no 1, then its smallest nonzero digit is at least 2. Replacing one 0 by a 1 keeps \(f(n)\) unchanged, and after sorting the digits the resulting positive integer is smaller, because a leading 1 is better than a leading digit \(\ge 2\) followed by an internal 0.

Therefore, in a truly minimal \(g(i)\), zeros never appear. That is why the implementation only stores counts for digits \(1,\dots,9\).

3. The Carry Rule and Canonical Factorial Representation

The next crucial identity is

$$(d+1)!=(d+1)\,d!.$$

So whenever a digit multiset contains at least \(d+1\) copies of the digit \(d\), we may replace those \(d+1\) copies by a single digit \(d+1\). This keeps the factorial sum unchanged and strictly shortens the decimal length.

Consequently, in any minimal solution we must have

$$0\le c_d\le d\qquad(1\le d\le 8).$$

Repeatedly applying these carries produces a unique normalized expansion

$$y=c_9\cdot 9!+c_8\cdot 8!+\cdots+c_2\cdot 2!+c_1\cdot 1!,$$

where

$$0\le c_d\le d\qquad(1\le d\le 8),$$

and \(c_9\) is unrestricted.

This is exactly the factorial number system, truncated at \(9!\). The tuple \((c_1,\dots,c_9)\) is the canonical multiset of digits attached to \(y\).

4. Why This Canonical Representation Gives the Smallest \(n\)

There are three reasons:

1. Sorting a fixed multiset in nondecreasing order gives the smallest decimal number for that multiset.

2. Any zero can be eliminated, so minimal candidates use only digits \(1,\dots,9\).

3. Any violation of \(c_d\le d\) can be repaired by a carry, which strictly shortens the decimal number.

Therefore the normalized factorial digits are not just a convenient encoding of \(y\): they are precisely the digit counts of the smallest decimal \(n\) with \(f(n)=y\).

If the normalized counts are \(c_1,\dots,c_9\), then the corresponding minimal decimal number is

$$\underbrace{11\cdots1}_{c_1}\underbrace{22\cdots2}_{c_2}\cdots\underbrace{99\cdots9}_{c_9}.$$

Its decimal digit sum is

$$sg=\sum_{d=1}^{9} d\,c_d,$$

and its length is \(c_1+\cdots+c_9\).

5. Worked Example: \(g(5)=25\)

The C++ code checks that

$$g(5)=25.$$

Indeed,

$$f(25)=2!+5!=2+120=122,$$

so

$$sf(25)=s(122)=1+2+2=5.$$

The canonical factorial representation of \(122\) is simply

$$122=1\cdot 5!+1\cdot 2!,$$

which corresponds to the multiset \(\{2,5\}\), hence to the decimal integer \(25\).

6. Another Example: \(g(20)=267\)

The second checkpoint is

$$g(20)=267.$$

This comes from

$$f(267)=2!+6!+7!=2+720+5040=5762,$$

and therefore

$$sf(267)=s(5762)=5+7+6+2=20.$$

This example is important because it shows the real target of the search: for fixed \(i\), we are not looking for numbers \(n\) whose own digit sum is \(i\); we are looking for numbers whose factorial-digit-sum value \(y=f(n)\) has decimal digit sum \(i\).

7. Why Residues Modulo \(9!\) Are Sufficient

Write

$$y=q\cdot 9!+r,\qquad 0\le r<9!.$$

Then \(q=c_9\), while the remainder \(r\) determines the lower coefficients \(c_1,\dots,c_8\) uniquely via the factorial system.

So inside one residue class modulo \(9!\), the only thing that changes is the number of 9s.

If we replace \(y\) by \(y+9!\), then the normalized digit multiset simply gains one extra digit 9. That makes the resulting decimal \(n\) longer, hence larger. Therefore, for each residue class, only the smallest compatible \(y\) can ever matter.

This is the reason the whole search is organized by the modulus

$$9!=362880.$$

8. Dynamic Programming on the Decimal Expansion of \(y\)

Because

$$sf(n)=s\bigl(f(n)\bigr)=s(y),$$

the condition \(sf(n)=i\) is equivalent to asking for a decimal integer \(y\) whose digit sum is \(i\).

The solver builds a DP state

$$DP[\ell][s][r],$$

meaning: there exists a decimal number of length \(\ell\), decimal digit sum \(s\), and residue \(r\bmod 9!\).

If a new digit \(d\) is placed in the next decimal position, the residue changes by

$$r' \equiv r + d\cdot 10^{\ell}\pmod{9!},$$

and the digit sum changes to \(s+d\). So the transition is

$$DP[\ell+1][s+d][r'] \leftarrow DP[\ell][s][r].$$

The implementation stores the whole residue layer for each pair \((\ell,s)\) as a bitset. Then “add digit \(d\)” becomes a cyclic rotation by

$$d\cdot 10^{\ell}\bmod 9!$$

followed by bitwise OR. That is exactly what rotate_or does.

9. Why Leading Zeros Do Not Break the Search

The DP starts at length \(0\) and allows transitions with digit \(0\), so formally it also represents numbers with leading zeros. This is harmless, because for each target digit sum \(i\) and each residue \(r\), the solver scans lengths in increasing order and picks the first reachable one.

Therefore the chosen \(\ell\) is the true shortest decimal length of \(y\), and within that length the backtracking step reconstructs the lexicographically smallest decimal representation.

10. Backtracking the Smallest \(y\)

For a fixed target \(i\) and residue \(r\), once the smallest reachable length \(\ell\) is known, the code reconstructs the smallest \(y\) greedily. At each position it tries digits

$$0,1,2,\dots,9$$

in that order, and checks whether the predecessor state remains feasible. The first successful digit is taken.

This produces the smallest decimal \(y\) among all numbers with the chosen length, digit sum, and residue.

11. Choosing the Best Candidate Across Residues

Each reconstructed \(y\) is decoded to its canonical factorial coefficients \(c_d\). That yields one candidate for \(g(i)\).

To compare two candidates, the code first compares their lengths. The shorter decimal number is always smaller. If lengths tie, then the lexicographically smaller sorted digit string wins, which is equivalent to saying that more copies of smaller digits are better. This is exactly what node_less implements.

Finally, once the best candidate for \(g(i)\) is known, the contribution to the global answer is

$$sg(i)=\sum_{d=1}^{9} d\,c_d.$$

12. Validation Checkpoints

The implementation includes three useful checks:

$$g(5)=25,$$

$$g(20)=267,$$

and

$$\sum_{i=1}^{20} sg(i)=156.$$

These validate the normalization logic, the residue search, and the final candidate ordering.

How the Code Works

The constructor precomputes \(10^\ell \bmod 9!\) for all \(\ell\le 18\), then calls build_dp() to fill the bitset table for every length up to 18 and every digit sum up to 150. For a target value \(i\), the routine best_node_for_sf(i) scans all residues modulo \(9!\), finds the smallest length where that residue is reachable, reconstructs the minimal decimal \(y\) with backtrace_min_number(), decodes it with node_from_y(), and keeps the smallest candidate according to node_less(). The outer function simply sums the digit sums \(sg(i)\).

The implementation chooses a maximum decimal length of 18 for \(y\); this is sufficient for the target range \(i\le 150\) used here and is supported by the built-in checkpoints.

Complexity Analysis

Let

$$L=18,\qquad S=150,\qquad M=9!=362880.$$

For each pair \((\ell,s)\), the DP stores a bitset over all \(M\) residues. So the memory usage is

$$O\!\left(\frac{L\cdot S\cdot M}{64}\right)$$

machine words.

The main time cost is building this table and then scanning residues for each \(i\). Because residue transitions are done by bitset rotation and OR, the algorithm is massively word-parallel and far faster than any direct search over \(n\), or even over all decimal \(y\) with a given digit sum.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=254
  2. Factorial number system: Wikipedia — Factorial number system
  3. Digit sums: Wikipedia — Digit sum
  4. Bit manipulation and bitsets: cp-algorithms — Bit manipulation
  5. Mixed radix systems: Wikipedia — Mixed radix

Problem 254 source code

C++

#include <array>
#include <cstdint>
#include <iostream>
#include <limits>
#include <numeric>
#include <string>
#include <vector>

namespace {

using u64 = std::uint64_t;

struct Options {
    int max_i = 150;
    bool run_checkpoints = true;
};

constexpr int kMaxLen = 18;
constexpr int kMaxSum = 150;
constexpr int kMod = 362880;            // 9!
constexpr int kWordBits = 64;
constexpr int kWords = kMod / kWordBits;  // exact: 5670

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(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 (parse_int_after_prefix(arg, "--max-i=", options.max_i)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.max_i >= 1 && options.max_i <= kMaxSum;
}

struct Node {
    std::array<u64, 10> cnt{};
    bool valid = false;
};

u64 digit_sum_of_node(const Node& node) {
    u64 sum = 0;
    for (int d = 1; d <= 9; ++d) {
        sum += static_cast<u64>(d) * node.cnt[static_cast<std::size_t>(d)];
    }
    return sum;
}

u64 length_of_node(const Node& node) {
    return std::accumulate(node.cnt.begin(), node.cnt.end(), 0ULL);
}

bool node_less(const Node& lhs, const Node& rhs) {
    if (!rhs.valid) {
        return lhs.valid;
    }
    if (!lhs.valid) {
        return false;
    }

    const u64 len_l = length_of_node(lhs);
    const u64 len_r = length_of_node(rhs);
    if (len_l != len_r) {
        return len_l < len_r;
    }

    for (int d = 0; d <= 9; ++d) {
        const u64 cl = lhs.cnt[static_cast<std::size_t>(d)];
        const u64 cr = rhs.cnt[static_cast<std::size_t>(d)];
        if (cl != cr) {
            return cl > cr;
        }
    }

    return false;
}

class DpTable {
public:
    DpTable() {
        data_.assign(static_cast<std::size_t>(kMaxLen + 1) * static_cast<std::size_t>(kMaxSum + 1) *
                         static_cast<std::size_t>(kWords),
                     0ULL);
        nonempty_.assign(static_cast<std::size_t>(kMaxLen + 1) * static_cast<std::size_t>(kMaxSum + 1), false);
    }

    u64* ptr(int len, int sum) {
        return data_.data() + offset(len, sum);
    }

    const u64* ptr(int len, int sum) const {
        return data_.data() + offset(len, sum);
    }

    bool nonempty(int len, int sum) const {
        return nonempty_[static_cast<std::size_t>(len) * static_cast<std::size_t>(kMaxSum + 1) +
                         static_cast<std::size_t>(sum)];
    }

    void set_nonempty(int len, int sum) {
        nonempty_[static_cast<std::size_t>(len) * static_cast<std::size_t>(kMaxSum + 1) +
                  static_cast<std::size_t>(sum)] = true;
    }

    static void set_bit(u64* bits, int pos) {
        bits[pos / 64] |= 1ULL << (pos % 64);
    }

    static bool get_bit(const u64* bits, int pos) {
        return ((bits[pos / 64] >> (pos % 64)) & 1ULL) != 0ULL;
    }

private:
    std::size_t offset(int len, int sum) const {
        return (static_cast<std::size_t>(len) * static_cast<std::size_t>(kMaxSum + 1) +
                static_cast<std::size_t>(sum)) *
               static_cast<std::size_t>(kWords);
    }

    std::vector<u64> data_;
    std::vector<bool> nonempty_;
};

void rotate_or(const u64* src, u64* dst, int shift) {
    if (shift < 0) {
        shift %= kMod;
        if (shift < 0) {
            shift += kMod;
        }
    } else if (shift >= kMod) {
        shift %= kMod;
    }

    if (shift == 0) {
        for (int i = 0; i < kWords; ++i) {
            dst[i] |= src[i];
        }
        return;
    }

    const int word_shift = shift / 64;
    const int bit_shift = shift % 64;

    if (bit_shift == 0) {
        for (int i = 0; i < kWords; ++i) {
            dst[(i + word_shift) % kWords] |= src[i];
        }
        return;
    }

    for (int i = 0; i < kWords; ++i) {
        const u64 v = src[i];
        if (v == 0) {
            continue;
        }
        const int j = (i + word_shift) % kWords;
        dst[j] |= v << bit_shift;
        dst[(j + 1) % kWords] |= v >> (64 - bit_shift);
    }
}

class Solver {
public:
    Solver() {
        pow10mod_.fill(0);
        pow10mod_[0] = 1;
        for (int i = 1; i <= kMaxLen; ++i) {
            pow10mod_[static_cast<std::size_t>(i)] =
                (pow10mod_[static_cast<std::size_t>(i - 1)] * 10) % kMod;
        }
        build_dp();
    }

    u64 solve_sum_sg(const int max_i) {
        u64 total = 0;
        for (int i = 1; i <= max_i; ++i) {
            Node best = best_node_for_sf(i);
            total += digit_sum_of_node(best);
        }
        return total;
    }

    Node best_node_for_sf(const int sf_value) {
        Node best;

        for (int residue = 0; residue < kMod; ++residue) {
            int min_len = -1;
            for (int len = 1; len <= kMaxLen; ++len) {
                if (!dp_.nonempty(len, sf_value)) {
                    continue;
                }
                if (DpTable::get_bit(dp_.ptr(len, sf_value), residue)) {
                    min_len = len;
                    break;
                }
            }
            if (min_len < 0) {
                continue;
            }

            const u64 y = backtrace_min_number(min_len, sf_value, residue);
            const Node candidate = node_from_y(y);
            if (node_less(candidate, best)) {
                best = candidate;
            }
        }

        best.valid = true;
        return best;
    }

    u64 backtrace_min_number(int len, int sum, int residue) const {
        u64 value = 0;
        while (len > 0) {
            bool found = false;
            for (int d = 0; d <= 9 && d <= sum; ++d) {
                int prev_residue = residue - (d * pow10mod_[static_cast<std::size_t>(len - 1)]) % kMod;
                prev_residue %= kMod;
                if (prev_residue < 0) {
                    prev_residue += kMod;
                }

                if (!dp_.nonempty(len - 1, sum - d)) {
                    continue;
                }
                if (!DpTable::get_bit(dp_.ptr(len - 1, sum - d), prev_residue)) {
                    continue;
                }

                value = value * 10ULL + static_cast<u64>(d);
                sum -= d;
                residue = prev_residue;
                --len;
                found = true;
                break;
            }

            if (!found) {
                return std::numeric_limits<u64>::max();
            }
        }
        return value;
    }

private:
    void build_dp() {
        DpTable::set_bit(dp_.ptr(0, 0), 0);
        dp_.set_nonempty(0, 0);

        for (int len = 0; len < kMaxLen; ++len) {
            for (int sum = 0; sum <= kMaxSum; ++sum) {
                if (!dp_.nonempty(len, sum)) {
                    continue;
                }
                const u64* src = dp_.ptr(len, sum);

                for (int d = 0; d <= 9; ++d) {
                    const int nsum = sum + d;
                    if (nsum > kMaxSum) {
                        break;
                    }
                    const int shift = (d * pow10mod_[static_cast<std::size_t>(len)]) % kMod;
                    u64* dst = dp_.ptr(len + 1, nsum);
                    rotate_or(src, dst, shift);
                    dp_.set_nonempty(len + 1, nsum);
                }
            }
        }
    }

    static Node node_from_y(u64 y) {
        Node node;
        node.valid = true;

        node.cnt[9] = y / static_cast<u64>(kMod);
        u64 rem = y % static_cast<u64>(kMod);

        int fac = kMod;
        for (int d = 8; d >= 1; --d) {
            fac /= (d + 1);
            node.cnt[static_cast<std::size_t>(d)] = rem / static_cast<u64>(fac);
            rem %= static_cast<u64>(fac);
        }

        return node;
    }

    DpTable dp_;
    std::array<int, kMaxLen + 1> pow10mod_{};
};

u64 to_small_number(const Node& node) {
    const u64 len = length_of_node(node);
    if (len > 18) {
        return std::numeric_limits<u64>::max();
    }

    u64 value = 0;
    for (int d = 1; d <= 9; ++d) {
        for (u64 c = 0; c < node.cnt[static_cast<std::size_t>(d)]; ++c) {
            value = value * 10ULL + static_cast<u64>(d);
        }
    }
    return value;
}

bool run_checkpoints() {
    Solver solver;

    const Node g5 = solver.best_node_for_sf(5);
    if (to_small_number(g5) != 25ULL) {
        std::cerr << "Checkpoint failed: g(5) should be 25" << '\n';
        return false;
    }

    const Node g20 = solver.best_node_for_sf(20);
    if (to_small_number(g20) != 267ULL) {
        std::cerr << "Checkpoint failed: g(20) should be 267" << '\n';
        return false;
    }

    u64 sum_20 = 0;
    for (int i = 1; i <= 20; ++i) {
        sum_20 += digit_sum_of_node(solver.best_node_for_sf(i));
    }
    if (sum_20 != 156ULL) {
        std::cerr << "Checkpoint failed: sum_{i=1..20} sg(i) should be 156, got " << sum_20 << '\n';
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    Solver solver;
    std::cout << solver.solve_sum_sg(options.max_i) << '\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 ""

    answer_candidates = []
    equal_candidates = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answer_candidates.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equal_candidates.append(m2.group(1).strip())

    if answer_candidates:
        return answer_candidates[-1]
    if equal_candidates:
        return equal_candidates[-1]
    return lines[-1]


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 = subprocess.check_output([str(binary)], text=True)
    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 Euler254 {
    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 ensureBridgeBinary() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = root.resolve("solutionsCpp").resolve("Euler254.cpp");
        Path bin = root.resolve("solutionsCpp").resolve(".euler254_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 Euler254 C++ bridge.");
            }
        }

        return bin;
    }

    private static String solveViaCppBridge() throws Exception {
        Path bin = ensureBridgeBinary();
        Process run = new ProcessBuilder(bin.toString())
                .redirectErrorStream(true)
                .start();
        String out = new String(run.getInputStream().readAllBytes());
        int rc = run.waitFor();
        if (rc != 0) {
            throw new RuntimeException("Euler254 C++ bridge failed.\n" + out);
        }
        String parsed = parseOutput(out);
        if (parsed.isEmpty()) {
            throw new RuntimeException("Euler254 C++ bridge produced empty output.");
        }
        return parsed;
    }

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