Problem 254: Sums of Digit Factorials

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$, let $f(n)$ be the sum of the factorials of its decimal digits, for example $f(342)=3!+4!+2!=32$, and let $sf(n)$ be the sum of the digits of $f(n)$, so $sf(342)=3+2=5$. Let $g(i)$ be the smallest positive integer $n$ with $sf(n)=i$, and let $sg(i)$ be the sum of the digits of $g(i)$. The statement notes that $sf(25)=5$ as well and that $g(5)=25$, so $sg(5)=2+5=7$; it also gives $g(20)=267$ and

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

We must find $\sum_{i=1}^{150} sg(i)$. ^[1] Below, $s(x)$ denotes the decimal digit sum of $x$, so that $sf(n)=s(f(n))$ and $sg(i)=s(g(i))$. The function $f$ is the OEIS sequence A061602; its fixed points $1$, $2$, $145$ and $40585$ are called factorions. ^[2]

The whole solution in one sentence: instead of trying the numbers $n$ one by one, we try the possible values $y$ of $f(n)$, because the smallest $n$ with $f(n)=y$ can be written down at once by dividing $y$ by $9!,8!,\dots,1!$ in turn, and only one value $y$ for each remainder after division by $9!=362880$ is worth trying, which a table of remainders finds quickly.

Mathematical Approach

We first compute small values and see why a search through $n$ fails. Sections 3 to 6 find the smallest $n$ with a given value of $f(n)$, Sections 7 and 8 leave one candidate for each remainder modulo $9!$, and Sections 9 and 10 find these candidates.

1. Small values by hand

The factorials of the digits are $0!=1!=1$, $2!=2$, $3!=6$, $4!=24$, $5!=120$, $6!=720$, $7!=5040$, $8!=40320$ and $9!=362880$. For the one-digit numbers $n=1,2,\dots,9$ the values of $sf(n)$ are therefore 1, 2, 6, 6, 3, 9, 9, 9 and 27; for instance $sf(5)=s(120)=3$ and $sf(9)=s(362880)=27$. This already gives $g(1)=1$, $g(2)=2$, $g(3)=5$, $g(6)=3$, $g(9)=6$ and $g(27)=9$. The value $5$ is missing, and the numbers $10,11,\dots,24$ give 2, 2, 3, 7, 7, 4, 10, 10, 10, 28, 3, 3, 4, 8 and 8. So $25$, with $f(25)=2!+5!=122$ and $s(122)=5$, is indeed the first number with $sf(n)=5$.

Continuing in the same way, the values $g(1),g(2),\dots,g(20)$ are 1, 2, 5, 15, 25, 3, 13, 23, 6, 16, 26, 44, 144, 256, 36, 136, 236, 67, 167 and 267. Their digit sums 1, 2, 5, 6, 7, 3, 4, 5, 6, 7, 8, 8, 9, 13, 9, 10, 11, 13, 14 and 15 add up to $156$, the value given in the statement.

2. Why a direct search cannot work

Scanning $n=1,2,3,\dots$ below $3\cdot10^6$ finds $g(i)$ exactly for $i\le44$; for instance $g(41)=2355679$ and $g(44)=2378889$, while $g(45)=12378889$ already has eight digits. Then the values explode: $g(63)$ has $42$ digits, $g(81)$ has $2779$, and $g(100)$ has $551165$ digits, $551146$ of them nines.

The reason is a lower bound. The smallest number with digit sum $i$ has as few digits as possible and, among those, the smallest leading digit, so all its other digits are nines. Writing $i=9u+v$ with $0\le v\le8$, it is the digit $v$, omitted when $v=0$, followed by $u$ nines:

$$y_0(i)=(v+1)\cdot10^{u}-1.$$

Every digit of $n$ contributes at most $9!$ to $f(n)$, and $f(g(i))\ge y_0(i)$ because $f(g(i))$ has digit sum $i$. Hence $g(i)$ has at least $\lceil y_0(i)/9!\rceil$ digits, which is $551147$ for $i=100$, since $y_0(100)/9!\approx551146.4$. We will never write $g(i)$ down: the task only needs its digit sum, which follows from the number of times each digit occurs.

3. Turning the search around

The value $sf(n)=s(f(n))$ depends on $n$ only through $y=f(n)$. For every positive integer $y$ let $m(y)$ be the smallest positive $n$ with $f(n)=y$; it exists, because $f$ maps the number written with $y$ ones to $y$. Then

$$g(i)=\min\{\,m(y) : s(y)=i\,\}.$$

Indeed, every $n$ with $sf(n)=i$ satisfies $n\ge m(f(n))$ and $s(f(n))=i$, and every $m(y)$ with $s(y)=i$ has $sf(m(y))=s(y)=i$. For example, $m(122)=25$ and $s(122)=5$. The work therefore splits into two questions: what is $m(y)$, and which $y$ with digit sum $i$ gives the smallest $m(y)$?

4. Digits as coins, and why zeros never help

If the digit $d$ occurs $c_d$ times in $n$, then $f(n)=\sum_d c_d\,d!$. Only the multiset of digits matters, as if each digit $d$ were a coin worth $d!$. A number with fewer digits is smaller, and among the arrangements of one multiset the smallest number lists the digits in increasing order, except that a zero cannot stand in front.

The smallest number $m(y)$ contains no zero. Since $0!=1!$, a zero is worth as much as a one, and there are three cases. Two zeros can be replaced by one $2$, because $1+1=2!$, and a zero together with a one can be replaced by a $2$ as well; both replacements keep $f$ and remove a digit. If there is exactly one zero and no one, the leading digit $a$ is at least $2$, and turning the zero into a one keeps $f$ and the length but lets the number start with $1\lt a$. For example, $f(102)=f(22)=4$ and $f(20)=f(12)=3$.

5. Carrying: the factorial number system

Because $(d+1)\cdot d!=(d+1)!$, any $d+1$ copies of the digit $d$ are worth exactly one digit $d+1$: $f(11)=f(2)$, $f(222)=f(3)$, $f(3333)=f(4)$, and so on. Replacing them removes $d$ digits, so in $m(y)$ each digit $d\le8$ occurs at most $d$ times. Only the nines are unlimited, because there is no digit $10$.

These two rules determine the digits completely. Suppose that $y=\sum_{d=1}^{9}c_d\,d!$ with $0\le c_d\le d$ for $d\le8$. Adding the identities $d\cdot d!=(d+1)!-d!$ gives the telescoping sum

$$\sum_{d=1}^{k}d\cdot d!=(k+1)!-1,$$

so the part $\sum_{d\le k}c_d\,d!$ is smaller than $(k+1)!$, while every term with $d\gt k$ is a multiple of $(k+1)!$. Hence that part is the remainder of $y$ modulo $(k+1)!$, and the counts are $c_9=\lfloor y/9!\rfloor$ and, for $1\le k\le8$,

$$c_k=\left\lfloor\frac{y\bmod (k+1)!}{k!}\right\rfloor.$$

This is the factorial number system, in which the digit at the place of $k!$ is at most $k$ and every natural number has exactly one representation; ^[3] here the top place $9!$ may hold any count, because decimal digits stop at $9$. Since $m(y)$ obeys both rules, we have proved the key fact.

Lemma. The number $m(y)$ consists of $c_1$ ones, $c_2$ twos, …, $c_9$ nines in increasing order, where $c_1,\dots,c_9$ are the factorial digits of $y$ above. It has $c_1+\dots+c_9$ digits, and its digit sum is $\sum_d d\,c_d$.

For example, $5762=7!+6!+2!$ gives $m(5762)=267$, and $362899=9!+3\cdot3!+1!$ gives $m(362899)=13339$. The counts are exactly what the greedy method of the change-making problem produces: take as many coins of the largest value as fit, then continue with the rest. The lemma says that for the coins $1!,2!,\dots,9!$ the greedy choice is the only one with the fewest coins, so this coin system is canonical. ^[4]

6. Comparing two candidates

By the lemma every candidate $m(y)$ is a sorted string of the digits $1$ to $9$. Of two such strings, the shorter one is the smaller number. If they have the same length, the one with more ones is smaller, because at the first position where they differ it has a $1$ and the other has a larger digit; if the numbers of ones agree, the one with more twos is smaller, and so on. So candidates are compared by length first and then by $c_1,c_2,\dots,c_9$ in turn, a larger count winning.

For $i=20$, no $m(y)$ with $s(y)=20$ has one or two digits, and exactly four have three digits: $267$, $269$, $278$ and $288$, from $y=5762$, $363602$, $45362$ and $80642$. All four have one $2$ and no $3$, $4$ or $5$; then $267$ and $269$ win with a $6$, and $267$ beats $269$ with a $7$. This is $g(20)=267$.

7. One candidate per remainder modulo 9!

Write $y=q\cdot9!+r$ with $0\le r\lt 9!$. By Section 5, $c_9=q$, while $c_1,\dots,c_8$ are the factorial digits of $r$ alone, since every $(k+1)!$ with $k\le8$ divides $9!$. So $m(y)$ has $q+L(r)$ digits, where $L(r)=c_1+\dots+c_8$ depends only on $r$ and is at most $1+2+\dots+8=36$. Within one class of remainders a larger $q$ means a longer number, so only the smallest $y\equiv r\pmod{9!}$ with $s(y)=i$ can give $g(i)$. This leaves at most one candidate for each of the $362880$ remainders $r$.

Far fewer remainders actually occur. A number and its digit sum leave the same remainder after division by $9$, because $10^j-1$ is divisible by $9$ for every $j$. ^[5] Since $9$ divides $9!$, every candidate has $r\equiv y\equiv i\pmod 9$, and only $9!/9=40320$ remainders are possible. For example, the class of $r=19$ consists of $19$, $362899$, $725779$, …; for $i=37$ its first member has digit sum $10$, but the second has digit sum $37$, so the candidate of this class is $m(362899)=13339$. Here $19\equiv37\equiv1\pmod 9$, as it must be.

8. Eighteen decimal digits are enough

For $i\le150$ there is a number $y\lt10^{17}$ with digit sum $i$, since $17$ digits can add up to $153$. By Section 7 its $m(y)$ has at most $\lfloor 10^{17}/9!\rfloor+36$ digits, that is $275573192275$. A value $y\ge10^{18}$ has $q\ge\lfloor10^{18}/9!\rfloor$, which is $2755731922398$, so its $m(y)$ is about ten times longer and can never be $g(i)$. Hence only values $y\lt10^{18}$, with at most $18$ decimal digits, need to be considered, and all of them fit into unsigned 64-bit integers, whose largest value is $2^{64}-1\approx1.8\cdot10^{19}$.

9. A table of reachable remainders

For each $i$ and each remainder $r$ we now need the smallest $y\equiv r\pmod{9!}$ with at most $18$ digits and $s(y)=i$. Dynamic programming solves such a task by building the answers to larger subproblems from stored answers to smaller ones. ^[6] Let $R(\ell,\sigma)$ be the set of remainders modulo $9!$ of the numbers below $10^{\ell}$, written as strings of $\ell$ digits with leading zeros allowed, whose digit sum is $\sigma$. Then $R(0,0)=\{0\}$, $R(0,\sigma)$ is empty for $\sigma\gt0$, and putting a new leading digit $d$ in front of a string of $\ell$ digits adds $d\cdot10^{\ell}$ to its value and $d$ to its digit sum:

$$\begin{aligned}&R(\ell+1,\sigma)\\ &\quad=\bigcup_{d=0}^{9}\Bigl(R(\ell,\sigma-d)+d\cdot10^{\ell}\Bigr),\end{aligned}$$

where adding a number to a set of remainders adds it to every element modulo $9!$, and the sets with $\sigma-d\lt0$ are empty. For instance, $R(1,\sigma)=\{\sigma\}$ for $\sigma\le9$, and $R(2,10)$ holds the nine remainders $19,28,37,\dots,91$, all congruent to $1$ modulo $9$. The smallest number of digits of a $y\equiv r$ with $s(y)=i$ is the smallest $\ell$ with $r\in R(\ell,i)$.

Each set is stored as a row of $9!$ bits, bit $r$ being $1$ when $r$ belongs to the set, so that a union is a bitwise OR of whole machine words. ^[7] Adding a constant $t$ to every element moves bit $r$ to bit $(r+t)\bmod 9!$, which is a cyclic rotation of the row. Since $9!=64\cdot5670$, a row is exactly $5670$ words of $64$ bits, and the rotation wraps around without a partly used word.

10. Rebuilding the smallest y digit by digit

Suppose $\ell$ is the smallest length with $r\in R(\ell,i)$. The smallest $y$ is then found from its leading digit downward: choose the smallest digit $d$ for which the rest can still be completed, that is, for which the remainder $(r-d\cdot10^{\ell-1})\bmod 9!$ belongs to $R(\ell-1,i-d)$, and continue in the same way with $\ell-1$ digits, digit sum $i-d$ and that remainder. Strings of equal length compare like numbers, by their first differing digit, so this greedy choice gives the smallest $y$; its leading digit is never $0$, because $\ell$ is minimal.

For $i=37$ and $r=19$ the smallest length is $\ell=6$: below $10^5$ the class of $19$ contains only $19$ itself, as $19+9!\gt10^5$, and $s(19)=10$. The leading digit $0$ would need $19\in R(5,37)$, which fails. The digits $1$ and $2$ would need the remainders $262899$ and $162899$, whose smallest members are already above $10^5$. The digit $3$ needs $62899\in R(5,34)$, and indeed $62899$ has digit sum $34=37-3$. Since $10^5\lt9!$, this five-digit rest is forced to be $62899$ itself, and the digits $6,2,8,9,9$ lead through the states $(\ell,\sigma,r)=(4,28,2899)$, $(3,26,899)$, $(2,18,99)$, $(1,9,9)$ and $(0,0,0)$. So $y=362899$.

Over all $40320$ admissible remainders for $i=37$, no candidate has fewer than five digits, and exactly four have five: $13339$, $13599$, $22479$ and $22499$, from $y=362899$, $725887$, $367948$ and $725788$. By Section 6, $13339$ wins: the last two have no $1$, and against $13599$, which also has one $1$ and no $2$, it has three $3$s instead of one. So $g(37)=13339$, in agreement with the direct search, and $sg(37)=19$.

How the Code Works

The Python and Java files are thin bridges to the C++ program. ^[8] Each compiles Euler254.cpp from the solutionsCpp folder with clang++ or g++ (-std=c++17 -O2) into a binary next to the source whenever that binary is missing or older than the source, runs it without arguments, so that the default limit and the checkpoints always apply, and returns the printed number. The parsers parse_output and parseOutput take the last value after answer: or an equals sign, else the last non-empty line, and reject empty output. Python finds the folder relative to its own file and Java relative to the working directory; a failed run raises an exception in Python and is detected through the exit code in Java. ^[9] All three tabs therefore run the computation described next.

The constants are kMaxLen = 18 (Section 8), kMaxSum = 150, kMod = 362880 and kWords = 5670. The class DpTable keeps all rows $R(\ell,\sigma)$ with $\ell\le18$ and $\sigma\le150$ in one flat vector of 64-bit words, with a flag nonempty per row for skipping empty rows. The constructor of Solver stores $10^{\ell}\bmod 9!$ in pow10mod_ and calls build_dp, which sets bit $0$ of row $(0,0)$ and then, for every nonempty row $(\ell,\sigma)$ with $\ell\lt18$ and every digit $d$ with $\sigma+d\le150$, ORs the row rotated by $d\cdot10^{\ell}\bmod 9!$ into row $(\ell+1,\sigma+d)$, as in Section 9. The helper rotate_or splits the shift into whole words and a bit offset below $64$ and ORs every nonzero source word into one or two neighbouring target words, taking the word index modulo kWords.

For a target value $i$, best_node_for_sf loops over all $9!$ remainders. For each one it finds the first length $\ell\le18$ whose row $(\ell,i)$ contains the remainder, and skips the remainder if there is none (Section 8). Otherwise backtrace_min_number rebuilds the smallest $y$ as in Section 10, and node_from_y turns $y$ into a Node of digit counts: cnt[9] is $\lfloor y/9!\rfloor$, and the factorial digits of the remainder fill cnt[8] down to cnt[1] (Section 5). node_less keeps the best candidate by the rule of Section 6, comparing length_of_node first and then the counts. Finally solve_sum_sg adds digit_sum_of_node, the sum $\sum_d d\,c_d$, over $i$ from $1$ to max_i; $g(i)$ itself is never written out.

In main, parse_arguments accepts --max-i=K with $1\le K\le150$, default 150, and --skip-checkpoints; other options and limits outside this range end the program with exit code 1, although a value of $K$ beyond the int range is not caught: its digits are collected in an int without an overflow check, which is undefined behaviour. Unless skipped, run_checkpoints builds its own Solver and checks $g(5)=25$, $g(20)=267$ and $\sum_{i\le20}sg(i)=156$, using to_small_number to write a node of at most 18 digits as an ordinary number; a failure ends the program with exit code 2. Then main builds a second Solver, so the table is built twice, and prints the sum. The values $y$ and the counts are uint64_t, which suffices by Section 8, and the remainders are ordinary int values below $9!$. The loop over remainders also visits the eight ninths of the remainders that violate $r\equiv i\pmod 9$ (Section 7); their bits are never set, so they cost time but change no result.

Complexity and Verification

The table has $19\cdot151\cdot5670=16\,267\,230$ words, about $130$ MB, which matches the memory the program uses. Only $1543$ of its $2869$ rows are ever filled, because $\ell$ digits have a digit sum of at most $9\ell$, and build_dp performs $13869$ rotations of $5670$ words each, about $7.9\cdot10^7$ word operations, which take a few hundredths of a second. For each $i$ the search makes at most $18$ bit tests for each of the $362880$ remainders and rebuilds at most $40320$ values $y$ with at most $18$ positions and $10$ digits per position: at most about $1.4\cdot10^7$ tests per value of $i$, or $2\cdot10^9$ for all $150$. The whole run, checkpoints included, takes about one and a half seconds. In general, with at most $D$ decimal digits, digit sums up to $S$, the modulus $M=9!$ and words of $w$ bits, the table needs $O(DSM)$ bits, building it takes $O(DSM/w)$ word operations, and the searches take $O(DSM)$ operations on single words.

Besides the three checkpoints of the program, the following checks were run for this article:

A brute force over all $n\lt3\cdot10^7$ gives $g(i)$ for every $i\le45$, including the values in Sections 1 and 2, in agreement with the lemma; its partial sum $\sum_{i\le45}sg(i)=659$ equals the output of the program run with --max-i=45.
A coin-change table confirms the lemma for every $y\lt726760$: the fewest digits from $1$ to $9$ with $f$-value $y$ number $c_1+\dots+c_9$, and only one multiset attains this. A brute force over $n\lt10^6$ confirms $m(y)$ for all $4014$ values $y$ whose $m(y)$ has at most six digits.
A search without the table of remainders: $m(y)$ has at least $\lfloor y/9!\rfloor$ digits and $m(y_0(i))$ at most $\lfloor y_0(i)/9!\rfloor+36$, so the $f$-value of $g(i)$ lies below $(\lfloor y_0(i)/9!\rfloor+37)\cdot9!$, and it suffices to decode every $y$ with digit sum $i$ from $y_0(i)$ up to that bound, about $1.3\cdot10^7$ numbers. For every $k=1,\dots,150$ its partial sum $\sum_{i\le k}sg(i)$ equals the output of the program run with --max-i=k, and its lengths of $g(63)$, $g(81)$ and $g(100)$ agree with the program.
A Python version of the table reproduces the $13869$ rotations and the $1543$ filled rows. Every reachable remainder satisfies $r\equiv i\pmod 9$; all $40320$ such remainders are reachable for $17\le i\le145$, but only $18$ for $i=1$, those of $1,10,\dots,10^{17}$.
The Python and Java bridges, run in a copy of the folder layout of the site, print the same number as the C++ program. This tests the bridges, not the method, since all three run the same C++ code.

Footnotes and References

The references below supply the background facts used above; the reduction to one candidate per remainder, the bounds and the operation counts are derived in this article.

Project Euler 254 — Sums of Digit Factorials. The official statement defines f(n), sf(n), g(i) and sg(i) as above, gives f(342) = 32, sf(342) = 5, g(5) = 25, sg(5) = 7, g(20) = 267 and the sum 156 of sg(i) for 1 ≤ i ≤ 20, and asks for the sum of sg(i) for 1 ≤ i ≤ 150.
OEIS A061602 — Sum of factorials of the digits of n. The entry lists f(n) from n = 0 on, with 0! = 1, and notes that the numbers with f(n) = n are called factorions and that in base 10 there are exactly four of them: 1, 2, 145 and 40585.
Wikipedia — Factorial number system. In this mixed-radix system the digit at the place of k! is at most k, and every natural number has exactly one representation; the article gives the identity that the sum of i·i! for i = 1, …, n equals (n + 1)! − 1 and finds the digits by repeated division with remainder.
Wikipedia — Change-making problem. The greedy method repeatedly takes the largest coin not greater than the remaining amount; a coin system is called canonical when this always gives the fewest coins, which fails for some systems, for example coins of 5, 10, 20 and 25 and the amount 40.
Wikipedia — Casting out nines. The digit sum of a number leaves the same remainder after division by 9 as the number itself, because the numbers 10ʲ − 1 are divisible by 9.
Wikipedia — Dynamic programming. The method needs optimal substructure and overlapping subproblems; the bottom-up approach solves the smaller subproblems first and builds on their stored solutions.
Wikipedia — Bit array. A bit array implements a set over {0, 1, …, n − 1}, a 1-bit marking an element that is present, and the union of two sets is computed word by word with OR in about n/w operations for words of w bits.
C++ — Euler254.cpp. The linked immutable C++ revision contains build_dp, rotate_or, best_node_for_sf, backtrace_min_number, node_from_y and the checks in run_checkpoints. It is the implementation record for this article.
Python — subprocess.check_output raises CalledProcessError when the called program returns a non-zero exit code; in Java, ProcessBuilder with redirectErrorStream(true) merges the standard error of the program into its standard output.

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());
    }
}