Problem 990: Addition Equations
View on Project EulerProject Euler Problem 990 Solution
EulerSolve provides an optimized solution for Project Euler Problem 990, Addition Equations, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Let \(F(n)\) denote the number of true decimal equations of the form \[ A_1+\cdots+A_\ell = B_1+\cdots+B_r, \] where every term is a positive integer written without leading zeros. The total length counts every digit, every plus sign, and the equality sign: \[ |E|=\sum_{i=1}^{\ell}\operatorname{len}(A_i)+\sum_{j=1}^{r}\operatorname{len}(B_j)+(\ell-1)+(r-1)+1. \] The problem asks for \(F(50)\bmod 10^9+7\). The implementations verify themselves against the small values \[ F(3)=9,\qquad F(5)=171,\qquad F(7)=4878. \] Brute force is not realistic: even before checking arithmetic, the number of possible strings explodes with the number of terms, their lengths, and all digit choices. The successful approach is to read the equation from right to left, one decimal column at a time, and count only the data that can still affect higher columns. Mathematical Approach The key simplification is that decimal addition is local. After the low columns have been fixed, the unfinished part of the equation is determined only by the terms that still have digits left and by the carry coming from below. Active Terms and the Length Bound A term is active in a column if it still contributes a digit in that column. If the equation has \(T=\ell+r\) total terms, then even the shortest possible equation has \(T\) digits and \(T-1\) separators, so \[ |E|\ge 2T-1....
Detailed mathematical approach
Problem Summary
Let \(F(n)\) denote the number of true decimal equations of the form
\[ A_1+\cdots+A_\ell = B_1+\cdots+B_r, \]
where every term is a positive integer written without leading zeros. The total length counts every digit, every plus sign, and the equality sign:
\[ |E|=\sum_{i=1}^{\ell}\operatorname{len}(A_i)+\sum_{j=1}^{r}\operatorname{len}(B_j)+(\ell-1)+(r-1)+1. \]
The problem asks for \(F(50)\bmod 10^9+7\). The implementations verify themselves against the small values
\[ F(3)=9,\qquad F(5)=171,\qquad F(7)=4878. \]
Brute force is not realistic: even before checking arithmetic, the number of possible strings explodes with the number of terms, their lengths, and all digit choices. The successful approach is to read the equation from right to left, one decimal column at a time, and count only the data that can still affect higher columns.
Mathematical Approach
The key simplification is that decimal addition is local. After the low columns have been fixed, the unfinished part of the equation is determined only by the terms that still have digits left and by the carry coming from below.
Active Terms and the Length Bound
A term is active in a column if it still contributes a digit in that column. If the equation has \(T=\ell+r\) total terms, then even the shortest possible equation has \(T\) digits and \(T-1\) separators, so
\[ |E|\ge 2T-1. \]
Therefore \(|E|\le 50\) implies \(T\le 25\). This is exactly the bound used by the implementations: there are never more than 25 active terms in total, and that keeps the state space finite and small.
Before any digit column is processed, the only committed characters are the separators. If there are \(\ell\) terms on the left and \(r\) on the right, the initial committed length is
\[ (\ell-1)+(r-1)+1=\ell+r-1. \]
One-Side Column Polynomials
Suppose a column begins with \(a\) active terms on one side, and exactly \(a'\le a\) of them survive to the next higher column. Then \(a-a'\) terms end in the current column, so their current digit is the leading digit of that number and must lie in \(\{1,\dots,9\}\). The surviving \(a'\) terms may contribute any digit in \(\{0,\dots,9\}\).
This gives the generating polynomial
\[ P_{a,a'}(x)=\bigl(x+x^2+\cdots+x^9\bigr)^{a-a'}\bigl(1+x+x^2+\cdots+x^9\bigr)^{a'}. \]
If we write
\[ c_{a,a'}(s)=[x^s]\,P_{a,a'}(x), \]
then \(c_{a,a'}(s)\) counts the admissible digit assignments on that side whose column sum is \(s\). Because the terms are distinguishable by their positions in the written equation, choosing which \(a'\) of the \(a\) active terms survive contributes the factor \(\binom{a}{a'}\).
The Carry Invariant
Let \(L\) and \(R\) be the left and right digit sums in the current column, and let \(\delta\) be the carry entering this column from the already processed lower columns. Valid base-10 arithmetic is equivalent to
\[ L-R+\delta=10\delta', \]
where \(\delta'\) is the carry passed to the next higher column. So only differences
\[ d=L-R \]
with \(d\equiv -\delta \pmod{10}\) are feasible for the current state.
For fixed survivor counts on both sides it is convenient to precompute
\[ T_{a,a',b,b'}(d)=\sum_{s-t=d} c_{a,a'}(s)\,c_{b,b'}(t), \]
the number of digit assignments that produce column difference \(d\). This is the real combinatorial object used by the implementations: once \(T_{a,a',b,b'}(d)\) is known, the individual digit sums \(s\) and \(t\) no longer have to be revisited during the DP.
State and Recurrence
Define \(D(u,a,b,\delta)\) as the number of partially constructed equations whose lower columns are already fixed, whose committed total length is \(u\), which still have \(a\) active left terms and \(b\) active right terms, and whose carry into the current column is \(\delta\).
The initial states are
\[ D(\ell+r-1,\ell,r,0)=1 \]
for every \(\ell\ge 1\), \(r\ge 1\), \(\ell+r\le 25\). Processing one new column appends one digit for every active term, so the length update is
\[ u'=u+a+b. \]
If the next column leaves \(a'\le a\) active terms on the left and \(b'\le b\) on the right, then the transition is
\[ D(u+a+b,a',b',\delta') \;+\!=\; D(u,a,b,\delta)\binom{a}{a'}\binom{b}{b'}T_{a,a',b,b'}(d), \]
whenever \(d+\delta=10\delta'\). An equation is complete exactly when
\[ a=b=0,\qquad \delta=0, \]
so the final count is
\[ F(n)=\sum_{u=0}^{n} D(u,0,0,0). \]
Why the Carry Window Is Only \([-25,25]\)
At any moment there are at most 25 active terms in total, so the column difference always satisfies
\[ |L-R|\le 9\cdot 25=225. \]
If \(|\delta|\le 25\), then
\[ |\delta'|=\left|\frac{L-R+\delta}{10}\right|\le \frac{225+25}{10}=25. \]
Since the process starts with carry \(0\), every reachable carry remains inside \([-25,25]\). That is why all three implementations can use a fixed carry dimension instead of an unbounded map.
Worked Example: Recovering \(F(5)=171\)
The length-5 check is a good sanity test because it already uses the same machinery as the full problem. Consider equations of the shape \(a+b=c\) with all three numbers one-digit. The separators contribute length 2, so one processed column gives total length 5. The initial state is \(D(2,2,1,0)=1\), and the only possible survivor counts are \(a'=0\), \(b'=0\).
The relevant polynomials are
\[ P_{2,0}(x)=\bigl(x+\cdots+x^9\bigr)^2,\qquad P_{1,0}(x)=x+\cdots+x^9. \]
Because the final carry must be zero, we need difference \(d=0\). For each right-hand digit \(c\in\{1,\dots,9\}\), there are exactly \(c-1\) ordered pairs \((a,b)\) with \(a+b=c\). Therefore
\[ T_{2,0,1,0}(0)=\sum_{c=2}^{9}(c-1)=36. \]
The mirrored shape \(c=a+b\) contributes another 36. Together with the 9 one-digit identities \(d=d\) and the 90 two-digit identities \(m=m\), we obtain
\[ F(5)=9+90+36+36=171, \]
which matches the checkpoint exactly.
How the Code Works
Precompute Combinatorial Tables
The C++, Python, and Java implementations first build all binomial coefficients up to 25. They then compute the coefficient lists \(c_{a,a'}(s)\) for every possible pair \((a,a')\) by multiplying the nonzero-leading-digit polynomial with the unrestricted-digit polynomial.
Build Difference Tables by Residue Class
For every quadruple \((a,a',b,b')\), the implementation combines the two one-side tables into the difference counts \(T_{a,a',b,b'}(d)\). These counts are grouped by \(d\bmod 10\), so a DP state with incoming carry \(\delta\) only has to inspect the compatible residue class \(d\equiv -\delta \pmod{10}\). This is the reason the transition step stays compact.
Sweep the DP and Read Off the Answer
The DP starts from every admissible split of the total terms between the two sides. States are processed in increasing committed length. Each update chooses survivor counts, multiplies by the corresponding binomial factors, reads the precomputed difference table, computes the next carry \((d+\delta)/10\), and adds the contribution to the state with length \(u+a+b\). After the sweep finishes, the answer is the sum of all completed states with length at most \(n\). The three implementations follow this same recurrence; the C++ version can additionally spread part of the preprocessing across multiple worker threads.
Complexity Analysis
The decisive fact is that \(|E|\le 50\) forces at most 25 total terms and only 51 possible carry values. The DP storage therefore has \((n+1)(26)(26)(51)\) slots; for \(n=50\) this is about \(1.76\times 10^6\) entries, and many of them are unreachable because \(a+b\le 25\).
The main work is the sweep over reachable states together with all survivor choices \((a',b')\) and the short precomputed lists of feasible differences in the matching residue class. Memory usage is linear in the number of stored states plus the finite family of difference tables. In practice this is fast because the algorithm never enumerates actual equation strings; it counts them indirectly through column sums, carries, and survivor sets.
Footnotes and References
- Problem page: Project Euler 990
- Positional notation: Wikipedia - Positional notation
- Carrying in arithmetic: Wikipedia - Carry (arithmetic)
- Generating functions: Wikipedia - Generating function
- Binomial coefficients: Wikipedia - Binomial coefficient
- Dynamic programming: Wikipedia - Dynamic programming
Problem 990 source code
C++
#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <string>
#include <thread>
#include <utility>
#include <vector>
namespace {
using u32 = std::uint32_t;
using u128 = unsigned __int128;
constexpr int MOD = 1'000'000'007;
constexpr int TARGET_N = 50;
constexpr int MAX_TERMS = (TARGET_N + 1) / 2;
constexpr int DELTA_LIMIT = 25;
constexpr int DELTA_SIZE = 2 * DELTA_LIMIT + 1;
struct Options {
bool run_checkpoints = true;
bool allow_multithreading = true;
unsigned requested_threads = 0;
int target_n = TARGET_N;
};
struct Transition {
std::array<std::vector<std::pair<int, u32>>, 10> by_residue;
};
int positive_mod(const int value, const int mod) {
const int result = value % mod;
return result < 0 ? result + mod : result;
}
u32 mod_add(const u32 a, const u32 b) {
const u32 sum = a + b;
return sum >= static_cast<u32>(MOD) ? sum - static_cast<u32>(MOD) : sum;
}
u32 mod_mul(const u32 a, const u32 b) {
return static_cast<u32>((static_cast<u128>(a) * static_cast<u128>(b)) % static_cast<u128>(MOD));
}
std::vector<u32> convolve(const std::vector<u32>& lhs, const std::vector<u32>& rhs) {
std::vector<u32> result(lhs.size() + rhs.size() - 1, 0);
for (std::size_t i = 0; i < lhs.size(); ++i) {
if (lhs[i] == 0) {
continue;
}
for (std::size_t j = 0; j < rhs.size(); ++j) {
if (rhs[j] == 0) {
continue;
}
const u32 add = mod_mul(lhs[i], rhs[j]);
result[i + j] = mod_add(result[i + j], add);
}
}
return result;
}
std::vector<u32> append_digit_set(const std::vector<u32>& poly, const int low_digit, const int high_digit) {
std::vector<u32> next(poly.size() + static_cast<std::size_t>(high_digit), 0);
for (std::size_t i = 0; i < poly.size(); ++i) {
if (poly[i] == 0) {
continue;
}
for (int digit = low_digit; digit <= high_digit; ++digit) {
next[i + static_cast<std::size_t>(digit)] =
mod_add(next[i + static_cast<std::size_t>(digit)], poly[i]);
}
}
return next;
}
struct Solver {
explicit Solver(const Options& options)
: thread_count(resolve_thread_count(options)),
binom(build_binom()),
digit_sums(build_digit_sum_polynomials()),
transition_index(build_transition_index()),
transitions(build_transitions()) {}
u32 solve(const int n) const {
assert(0 <= n && n <= TARGET_N);
const int stride_len = (MAX_TERMS + 1) * (MAX_TERMS + 1) * DELTA_SIZE;
const int stride_a = (MAX_TERMS + 1) * DELTA_SIZE;
const int stride_b = DELTA_SIZE;
const auto index = [&](const int used, const int left_active, const int right_active, const int delta) {
return used * stride_len + left_active * stride_a + right_active * stride_b +
(delta + DELTA_LIMIT);
};
std::vector<u32> dp(static_cast<std::size_t>(n + 1) * static_cast<std::size_t>(stride_len), 0);
for (int left_terms = 1; left_terms <= MAX_TERMS; ++left_terms) {
for (int right_terms = 1; right_terms + left_terms <= MAX_TERMS; ++right_terms) {
const int base_length = left_terms + right_terms - 1;
if (base_length > n) {
continue;
}
dp[static_cast<std::size_t>(index(base_length, left_terms, right_terms, 0))] = 1;
}
}
for (int used = 0; used <= n; ++used) {
for (int left_active = 0; left_active <= MAX_TERMS; ++left_active) {
for (int right_active = 0; right_active + left_active <= MAX_TERMS; ++right_active) {
if (left_active == 0 && right_active == 0) {
continue;
}
const int next_length = used + left_active + right_active;
if (next_length > n) {
continue;
}
for (int delta = -DELTA_LIMIT; delta <= DELTA_LIMIT; ++delta) {
const u32 current =
dp[static_cast<std::size_t>(index(used, left_active, right_active, delta))];
if (current == 0) {
continue;
}
for (int left_next = 0; left_next <= left_active; ++left_next) {
const u32 choose_left = binom[left_active][left_next];
for (int right_next = 0; right_next <= right_active; ++right_next) {
const u32 choose_right = binom[right_active][right_next];
const u32 structure_weight =
mod_mul(current, mod_mul(choose_left, choose_right));
const int tindex =
transition_index[left_active][left_next][right_active][right_next];
const auto& group =
transitions[static_cast<std::size_t>(tindex)]
.by_residue[static_cast<std::size_t>(positive_mod(-delta, 10))];
for (const auto& [diff, ways] : group) {
const int total = diff + delta;
const int next_delta = total / 10;
if (next_delta < -DELTA_LIMIT || next_delta > DELTA_LIMIT) {
continue;
}
const std::size_t target =
static_cast<std::size_t>(index(next_length,
left_next,
right_next,
next_delta));
dp[target] = mod_add(dp[target], mod_mul(structure_weight, ways));
}
}
}
}
}
}
}
u32 answer = 0;
for (int used = 0; used <= n; ++used) {
answer = mod_add(answer, dp[static_cast<std::size_t>(index(used, 0, 0, 0))]);
}
return answer;
}
private:
struct Task {
int left_active = 0;
int left_next = 0;
int right_active = 0;
int right_next = 0;
int index = -1;
};
unsigned thread_count;
std::array<std::array<u32, MAX_TERMS + 1>, MAX_TERMS + 1> binom{};
std::array<std::array<std::vector<u32>, MAX_TERMS + 1>, MAX_TERMS + 1> digit_sums{};
std::array<std::array<std::array<std::array<int, MAX_TERMS + 1>, MAX_TERMS + 1>, MAX_TERMS + 1>,
MAX_TERMS + 1>
transition_index{};
std::vector<Transition> transitions;
static unsigned resolve_thread_count(const Options& options) {
if (!options.allow_multithreading) {
return 1;
}
if (options.requested_threads != 0) {
return std::max(1U, options.requested_threads);
}
const unsigned detected = std::thread::hardware_concurrency();
return detected == 0 ? 1U : detected;
}
static std::array<std::array<u32, MAX_TERMS + 1>, MAX_TERMS + 1> build_binom() {
std::array<std::array<u32, MAX_TERMS + 1>, MAX_TERMS + 1> result{};
result[0][0] = 1;
for (int n = 1; n <= MAX_TERMS; ++n) {
result[n][0] = 1;
result[n][n] = 1;
for (int k = 1; k < n; ++k) {
result[n][k] = result[n - 1][k - 1] + result[n - 1][k];
}
}
return result;
}
static std::array<std::array<std::vector<u32>, MAX_TERMS + 1>, MAX_TERMS + 1>
build_digit_sum_polynomials() {
std::array<std::array<std::vector<u32>, MAX_TERMS + 1>, MAX_TERMS + 1> result{};
std::array<std::vector<u32>, MAX_TERMS + 1> free_digits{};
std::array<std::vector<u32>, MAX_TERMS + 1> leading_digits{};
free_digits[0] = {1};
leading_digits[0] = {1};
for (int terms = 1; terms <= MAX_TERMS; ++terms) {
free_digits[terms] = append_digit_set(free_digits[terms - 1], 0, 9);
leading_digits[terms] = append_digit_set(leading_digits[terms - 1], 1, 9);
}
for (int active = 0; active <= MAX_TERMS; ++active) {
for (int next = 0; next <= active; ++next) {
result[active][next] = convolve(leading_digits[active - next], free_digits[next]);
}
}
return result;
}
std::array<std::array<std::array<std::array<int, MAX_TERMS + 1>, MAX_TERMS + 1>, MAX_TERMS + 1>,
MAX_TERMS + 1>
build_transition_index() const {
std::array<std::array<std::array<std::array<int, MAX_TERMS + 1>, MAX_TERMS + 1>, MAX_TERMS + 1>,
MAX_TERMS + 1>
result{};
int next_index = 0;
for (int left_active = 0; left_active <= MAX_TERMS; ++left_active) {
for (int left_next = 0; left_next <= left_active; ++left_next) {
for (int right_active = 0; right_active + left_active <= MAX_TERMS; ++right_active) {
for (int right_next = 0; right_next <= right_active; ++right_next) {
result[left_active][left_next][right_active][right_next] = next_index++;
}
}
}
}
return result;
}
Transition build_single_transition(const Task& task) const {
Transition transition;
const auto& left_poly = digit_sums[task.left_active][task.left_next];
const auto& right_poly = digit_sums[task.right_active][task.right_next];
const int min_diff = -static_cast<int>(right_poly.size()) + 1;
const int max_diff = static_cast<int>(left_poly.size()) - 1;
std::vector<u32> diff_counts(static_cast<std::size_t>(max_diff - min_diff + 1), 0);
for (int left_sum = 0; left_sum < static_cast<int>(left_poly.size()); ++left_sum) {
if (left_poly[static_cast<std::size_t>(left_sum)] == 0) {
continue;
}
for (int right_sum = 0; right_sum < static_cast<int>(right_poly.size()); ++right_sum) {
if (right_poly[static_cast<std::size_t>(right_sum)] == 0) {
continue;
}
const int diff = left_sum - right_sum;
const u32 ways =
mod_mul(left_poly[static_cast<std::size_t>(left_sum)],
right_poly[static_cast<std::size_t>(right_sum)]);
diff_counts[static_cast<std::size_t>(diff - min_diff)] =
mod_add(diff_counts[static_cast<std::size_t>(diff - min_diff)], ways);
}
}
for (int diff = min_diff; diff <= max_diff; ++diff) {
const u32 ways = diff_counts[static_cast<std::size_t>(diff - min_diff)];
if (ways == 0) {
continue;
}
transition.by_residue[static_cast<std::size_t>(positive_mod(diff, 10))].push_back({diff, ways});
}
return transition;
}
std::vector<Transition> build_transitions() const {
std::vector<Task> tasks;
for (int left_active = 0; left_active <= MAX_TERMS; ++left_active) {
for (int left_next = 0; left_next <= left_active; ++left_next) {
for (int right_active = 0; right_active + left_active <= MAX_TERMS; ++right_active) {
for (int right_next = 0; right_next <= right_active; ++right_next) {
tasks.push_back({left_active,
left_next,
right_active,
right_next,
transition_index[left_active][left_next][right_active][right_next]});
}
}
}
}
std::vector<Transition> result(tasks.size());
const unsigned workers =
std::min<unsigned>(thread_count, static_cast<unsigned>(std::max<std::size_t>(1, tasks.size())));
if (workers <= 1) {
for (const Task& task : tasks) {
result[static_cast<std::size_t>(task.index)] = build_single_transition(task);
}
return result;
}
std::vector<std::thread> pool;
pool.reserve(workers);
for (unsigned worker = 0; worker < workers; ++worker) {
pool.emplace_back([&, worker]() {
for (std::size_t i = worker; i < tasks.size(); i += workers) {
const Task& task = tasks[i];
result[static_cast<std::size_t>(task.index)] = build_single_transition(task);
}
});
}
for (std::thread& thread : pool) {
thread.join();
}
return result;
}
};
void usage() {
std::cerr
<< "Usage:\n"
<< " ./Euler990 [n] [--skip-checkpoints] [--single-thread] [--threads=N]\n";
}
Options parse_options(const 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 (arg == "--single-thread") {
options.allow_multithreading = false;
options.requested_threads = 1;
continue;
}
if (arg.rfind("--threads=", 0) == 0) {
options.requested_threads = static_cast<unsigned>(std::stoul(arg.substr(10)));
continue;
}
if (!arg.empty() && arg[0] == '-') {
usage();
std::exit(EXIT_FAILURE);
}
options.target_n = std::stoi(arg);
}
if (options.target_n < 0 || options.target_n > TARGET_N) {
std::cerr << "n must satisfy 0 <= n <= " << TARGET_N << ".\n";
std::exit(EXIT_FAILURE);
}
return options;
}
void run_checkpoints(const Solver& solver) {
assert(solver.solve(0) == 0);
assert(solver.solve(1) == 0);
assert(solver.solve(2) == 0);
assert(solver.solve(3) == 9);
assert(solver.solve(5) == 171);
assert(solver.solve(7) == 4878);
}
} // namespace
int main(int argc, char** argv) {
const Options options = parse_options(argc, argv);
const Solver solver(options);
if (options.run_checkpoints) {
run_checkpoints(solver);
}
std::cout << solver.solve(options.target_n) << '\n';
return 0;
}
Python
from __future__ import annotations
import sys
from array import array
MOD = 1_000_000_007
TARGET_N = 50
MAX_TERMS = (TARGET_N + 1) // 2
DELTA_LIMIT = 25
DELTA_SIZE = 2 * DELTA_LIMIT + 1
class Options:
__slots__ = ("run_checkpoints", "allow_multithreading", "requested_threads", "target_n")
def __init__(self) -> None:
self.run_checkpoints = True
self.allow_multithreading = True
self.requested_threads = 0
self.target_n = TARGET_N
class Transition:
__slots__ = ("by_residue",)
def __init__(self) -> None:
self.by_residue = [[] for _ in range(10)]
def mod_add(a: int, b: int) -> int:
s = a + b
return s - MOD if s >= MOD else s
def mod_mul(a: int, b: int) -> int:
return (a * b) % MOD
def convolve(lhs: list[int], rhs: list[int]) -> list[int]:
result = [0] * (len(lhs) + len(rhs) - 1)
for i, left in enumerate(lhs):
if left == 0:
continue
for j, right in enumerate(rhs):
if right == 0:
continue
result[i + j] = mod_add(result[i + j], mod_mul(left, right))
return result
def append_digit_set(poly: list[int], low_digit: int, high_digit: int) -> list[int]:
next_poly = [0] * (len(poly) + high_digit)
for i, coeff in enumerate(poly):
if coeff == 0:
continue
for digit in range(low_digit, high_digit + 1):
next_poly[i + digit] = mod_add(next_poly[i + digit], coeff)
return next_poly
class Solver:
def __init__(self, _options: Options) -> None:
self.binom = self.build_binom()
self.digit_sums = self.build_digit_sum_polynomials()
self.transition_base, self.transition_count = self.build_transition_base()
self.transitions = self.build_transitions()
@staticmethod
def build_binom() -> list[list[int]]:
result = [[0] * (MAX_TERMS + 1) for _ in range(MAX_TERMS + 1)]
result[0][0] = 1
for n in range(1, MAX_TERMS + 1):
result[n][0] = 1
result[n][n] = 1
for k in range(1, n):
result[n][k] = result[n - 1][k - 1] + result[n - 1][k]
return result
@staticmethod
def build_digit_sum_polynomials() -> list[list[list[int] | None]]:
result: list[list[list[int] | None]] = [[None] * (MAX_TERMS + 1) for _ in range(MAX_TERMS + 1)]
free_digits: list[list[int] | None] = [None] * (MAX_TERMS + 1)
leading_digits: list[list[int] | None] = [None] * (MAX_TERMS + 1)
free_digits[0] = [1]
leading_digits[0] = [1]
for terms in range(1, MAX_TERMS + 1):
free_digits[terms] = append_digit_set(free_digits[terms - 1], 0, 9)
leading_digits[terms] = append_digit_set(leading_digits[terms - 1], 1, 9)
for active in range(MAX_TERMS + 1):
for nxt in range(active + 1):
result[active][nxt] = convolve(leading_digits[active - nxt], free_digits[nxt])
return result
@staticmethod
def build_transition_base() -> tuple[list[list[list[int]]], int]:
base = [[[0] * (MAX_TERMS + 1) for _ in range(MAX_TERMS + 1)] for _ in range(MAX_TERMS + 1)]
next_index = 0
for left_active in range(MAX_TERMS + 1):
for left_next in range(left_active + 1):
for right_active in range(MAX_TERMS - left_active + 1):
base[left_active][left_next][right_active] = next_index
next_index += right_active + 1
return base, next_index
def build_single_transition(
self,
left_active: int,
left_next: int,
right_active: int,
right_next: int,
) -> Transition:
transition = Transition()
left_poly = self.digit_sums[left_active][left_next]
right_poly = self.digit_sums[right_active][right_next]
min_diff = -len(right_poly) + 1
max_diff = len(left_poly) - 1
diff_counts = [0] * (max_diff - min_diff + 1)
for left_sum, left_ways in enumerate(left_poly):
if left_ways == 0:
continue
for right_sum, right_ways in enumerate(right_poly):
if right_ways == 0:
continue
diff = left_sum - right_sum
diff_counts[diff - min_diff] = mod_add(
diff_counts[diff - min_diff],
mod_mul(left_ways, right_ways),
)
for diff in range(min_diff, max_diff + 1):
ways = diff_counts[diff - min_diff]
if ways == 0:
continue
transition.by_residue[diff % 10].append((diff, ways))
return transition
def build_transitions(self) -> list[Transition]:
transitions: list[Transition | None] = [None] * self.transition_count
for left_active in range(MAX_TERMS + 1):
for left_next in range(left_active + 1):
for right_active in range(MAX_TERMS - left_active + 1):
base = self.transition_base[left_active][left_next][right_active]
for right_next in range(right_active + 1):
transitions[base + right_next] = self.build_single_transition(
left_active, left_next, right_active, right_next
)
return transitions
def solve(self, n: int) -> int:
assert 0 <= n <= TARGET_N
stride_len = (MAX_TERMS + 1) * (MAX_TERMS + 1) * DELTA_SIZE
stride_a = (MAX_TERMS + 1) * DELTA_SIZE
stride_b = DELTA_SIZE
def index(used: int, left_active: int, right_active: int, delta: int) -> int:
return (
used * stride_len
+ left_active * stride_a
+ right_active * stride_b
+ delta
+ DELTA_LIMIT
)
dp = array("I", [0]) * ((n + 1) * stride_len)
for left_terms in range(1, MAX_TERMS + 1):
max_right_terms = MAX_TERMS - left_terms
for right_terms in range(1, max_right_terms + 1):
base_length = left_terms + right_terms - 1
if base_length > n:
continue
dp[index(base_length, left_terms, right_terms, 0)] = 1
binom = self.binom
transition_base = self.transition_base
transitions = self.transitions
for used in range(n + 1):
for left_active in range(MAX_TERMS + 1):
max_right_active = MAX_TERMS - left_active
for right_active in range(max_right_active + 1):
if left_active == 0 and right_active == 0:
continue
next_length = used + left_active + right_active
if next_length > n:
continue
for delta in range(-DELTA_LIMIT, DELTA_LIMIT + 1):
current = dp[index(used, left_active, right_active, delta)]
if current == 0:
continue
residue = (-delta) % 10
for left_next in range(left_active + 1):
choose_left = binom[left_active][left_next]
for right_next in range(right_active + 1):
structure_weight = mod_mul(
current,
mod_mul(choose_left, binom[right_active][right_next]),
)
tindex = transition_base[left_active][left_next][right_active] + right_next
group = transitions[tindex].by_residue[residue]
for diff, ways in group:
total = diff + delta
next_delta = total // 10
if -DELTA_LIMIT <= next_delta <= DELTA_LIMIT:
target = index(next_length, left_next, right_next, next_delta)
dp[target] = mod_add(dp[target], mod_mul(structure_weight, ways))
answer = 0
for used in range(n + 1):
answer = mod_add(answer, dp[index(used, 0, 0, 0)])
return answer
def usage() -> None:
sys.stderr.write(
"Usage:\n"
" python Euler990.py [n] [--skip-checkpoints] [--single-thread] [--threads=N]\n"
)
def parse_options(argv: list[str]) -> Options:
options = Options()
for arg in argv[1:]:
if arg == "--skip-checkpoints":
options.run_checkpoints = False
continue
if arg == "--single-thread":
options.allow_multithreading = False
options.requested_threads = 1
continue
if arg.startswith("--threads="):
options.requested_threads = int(arg[10:])
continue
if arg.startswith("-"):
usage()
raise SystemExit(1)
options.target_n = int(arg)
if options.target_n < 0 or options.target_n > TARGET_N:
sys.stderr.write(f"n must satisfy 0 <= n <= {TARGET_N}.\n")
raise SystemExit(1)
return options
def run_checkpoints(solver: Solver) -> None:
assert solver.solve(0) == 0
assert solver.solve(1) == 0
assert solver.solve(2) == 0
assert solver.solve(3) == 9
assert solver.solve(5) == 171
assert solver.solve(7) == 4878
def main(argv: list[str]) -> int:
options = parse_options(argv)
solver = Solver(options)
if options.run_checkpoints:
run_checkpoints(solver)
print(solver.solve(options.target_n))
return 0
if __name__ == "__main__":
raise SystemExit(main(sys.argv))
Java
public class Euler990 {
private static final int MOD = 1_000_000_007;
private static final int TARGET_N = 50;
private static final int MAX_TERMS = (TARGET_N + 1) / 2;
private static final int DELTA_LIMIT = 25;
private static final int DELTA_SIZE = 2 * DELTA_LIMIT + 1;
private static final class Options {
boolean runCheckpoints = true;
boolean allowMultithreading = true;
int requestedThreads = 0;
int targetN = TARGET_N;
}
private static final class IntList {
private int[] data = new int[16];
private int size = 0;
void add(int value) {
if (size == data.length) {
int[] next = new int[data.length * 2];
System.arraycopy(data, 0, next, 0, data.length);
data = next;
}
data[size++] = value;
}
int[] toArray() {
int[] result = new int[size];
System.arraycopy(data, 0, result, 0, size);
return result;
}
}
private static final class Transition {
final int[][] diffs = new int[10][];
final int[][] ways = new int[10][];
}
private static int modAdd(int a, int b) {
int s = a + b;
return s >= MOD ? s - MOD : s;
}
private static int modMul(int a, int b) {
return (int) (((long) a * (long) b) % MOD);
}
private static int[] convolve(int[] lhs, int[] rhs) {
int[] result = new int[lhs.length + rhs.length - 1];
for (int i = 0; i < lhs.length; ++i) {
int left = lhs[i];
if (left == 0) {
continue;
}
for (int j = 0; j < rhs.length; ++j) {
int right = rhs[j];
if (right == 0) {
continue;
}
result[i + j] = modAdd(result[i + j], modMul(left, right));
}
}
return result;
}
private static int[] appendDigitSet(int[] poly, int lowDigit, int highDigit) {
int[] next = new int[poly.length + highDigit];
for (int i = 0; i < poly.length; ++i) {
int coeff = poly[i];
if (coeff == 0) {
continue;
}
for (int digit = lowDigit; digit <= highDigit; ++digit) {
next[i + digit] = modAdd(next[i + digit], coeff);
}
}
return next;
}
private static final class Solver {
private final int[][] binom;
private final int[][][] digitSums;
private final int[][][] transitionBase;
private final int transitionCount;
private final Transition[] transitions;
Solver(Options options) {
binom = buildBinom();
digitSums = buildDigitSumPolynomials();
transitionBase = new int[MAX_TERMS + 1][MAX_TERMS + 1][MAX_TERMS + 1];
transitionCount = buildTransitionBase();
transitions = buildTransitions();
}
private static int[][] buildBinom() {
int[][] result = new int[MAX_TERMS + 1][MAX_TERMS + 1];
result[0][0] = 1;
for (int n = 1; n <= MAX_TERMS; ++n) {
result[n][0] = 1;
result[n][n] = 1;
for (int k = 1; k < n; ++k) {
result[n][k] = result[n - 1][k - 1] + result[n - 1][k];
}
}
return result;
}
private static int[][][] buildDigitSumPolynomials() {
int[][][] result = new int[MAX_TERMS + 1][MAX_TERMS + 1][];
int[][] freeDigits = new int[MAX_TERMS + 1][];
int[][] leadingDigits = new int[MAX_TERMS + 1][];
freeDigits[0] = new int[] {1};
leadingDigits[0] = new int[] {1};
for (int terms = 1; terms <= MAX_TERMS; ++terms) {
freeDigits[terms] = appendDigitSet(freeDigits[terms - 1], 0, 9);
leadingDigits[terms] = appendDigitSet(leadingDigits[terms - 1], 1, 9);
}
for (int active = 0; active <= MAX_TERMS; ++active) {
for (int next = 0; next <= active; ++next) {
result[active][next] = convolve(leadingDigits[active - next], freeDigits[next]);
}
}
return result;
}
private int buildTransitionBase() {
int nextIndex = 0;
for (int leftActive = 0; leftActive <= MAX_TERMS; ++leftActive) {
for (int leftNext = 0; leftNext <= leftActive; ++leftNext) {
for (int rightActive = 0; rightActive + leftActive <= MAX_TERMS; ++rightActive) {
transitionBase[leftActive][leftNext][rightActive] = nextIndex;
nextIndex += rightActive + 1;
}
}
}
return nextIndex;
}
private Transition buildSingleTransition(int leftActive, int leftNext, int rightActive, int rightNext) {
Transition transition = new Transition();
int[] leftPoly = digitSums[leftActive][leftNext];
int[] rightPoly = digitSums[rightActive][rightNext];
int minDiff = -rightPoly.length + 1;
int maxDiff = leftPoly.length - 1;
int[] diffCounts = new int[maxDiff - minDiff + 1];
for (int leftSum = 0; leftSum < leftPoly.length; ++leftSum) {
int leftWays = leftPoly[leftSum];
if (leftWays == 0) {
continue;
}
for (int rightSum = 0; rightSum < rightPoly.length; ++rightSum) {
int rightWays = rightPoly[rightSum];
if (rightWays == 0) {
continue;
}
int diff = leftSum - rightSum;
diffCounts[diff - minDiff] =
modAdd(diffCounts[diff - minDiff], modMul(leftWays, rightWays));
}
}
IntList[] diffLists = new IntList[10];
IntList[] wayLists = new IntList[10];
for (int diff = minDiff; diff <= maxDiff; ++diff) {
int ways = diffCounts[diff - minDiff];
if (ways == 0) {
continue;
}
int residue = ((diff % 10) + 10) % 10;
if (diffLists[residue] == null) {
diffLists[residue] = new IntList();
wayLists[residue] = new IntList();
}
diffLists[residue].add(diff);
wayLists[residue].add(ways);
}
for (int residue = 0; residue < 10; ++residue) {
transition.diffs[residue] = diffLists[residue] == null ? new int[0] : diffLists[residue].toArray();
transition.ways[residue] = wayLists[residue] == null ? new int[0] : wayLists[residue].toArray();
}
return transition;
}
private Transition[] buildTransitions() {
Transition[] result = new Transition[transitionCount];
for (int leftActive = 0; leftActive <= MAX_TERMS; ++leftActive) {
for (int leftNext = 0; leftNext <= leftActive; ++leftNext) {
for (int rightActive = 0; rightActive + leftActive <= MAX_TERMS; ++rightActive) {
int base = transitionBase[leftActive][leftNext][rightActive];
for (int rightNext = 0; rightNext <= rightActive; ++rightNext) {
result[base + rightNext] =
buildSingleTransition(leftActive, leftNext, rightActive, rightNext);
}
}
}
}
return result;
}
int solve(int n) {
assert 0 <= n && n <= TARGET_N;
final int strideLen = (MAX_TERMS + 1) * (MAX_TERMS + 1) * DELTA_SIZE;
final int strideA = (MAX_TERMS + 1) * DELTA_SIZE;
final int strideB = DELTA_SIZE;
int[] dp = new int[(n + 1) * strideLen];
for (int leftTerms = 1; leftTerms <= MAX_TERMS; ++leftTerms) {
for (int rightTerms = 1; rightTerms + leftTerms <= MAX_TERMS; ++rightTerms) {
int baseLength = leftTerms + rightTerms - 1;
if (baseLength > n) {
continue;
}
dp[index(baseLength, leftTerms, rightTerms, 0, strideLen, strideA, strideB)] = 1;
}
}
for (int used = 0; used <= n; ++used) {
for (int leftActive = 0; leftActive <= MAX_TERMS; ++leftActive) {
for (int rightActive = 0; rightActive + leftActive <= MAX_TERMS; ++rightActive) {
if (leftActive == 0 && rightActive == 0) {
continue;
}
int nextLength = used + leftActive + rightActive;
if (nextLength > n) {
continue;
}
for (int delta = -DELTA_LIMIT; delta <= DELTA_LIMIT; ++delta) {
int current =
dp[index(used, leftActive, rightActive, delta, strideLen, strideA, strideB)];
if (current == 0) {
continue;
}
int residue = ((-delta) % 10 + 10) % 10;
for (int leftNext = 0; leftNext <= leftActive; ++leftNext) {
int chooseLeft = binom[leftActive][leftNext];
for (int rightNext = 0; rightNext <= rightActive; ++rightNext) {
int structureWeight =
modMul(current, modMul(chooseLeft, binom[rightActive][rightNext]));
int tindex = transitionBase[leftActive][leftNext][rightActive] + rightNext;
Transition transition = transitions[tindex];
int[] diffs = transition.diffs[residue];
int[] ways = transition.ways[residue];
for (int i = 0; i < diffs.length; ++i) {
int total = diffs[i] + delta;
int nextDelta = total / 10;
if (nextDelta < -DELTA_LIMIT || nextDelta > DELTA_LIMIT) {
continue;
}
int target =
index(nextLength, leftNext, rightNext, nextDelta, strideLen, strideA,
strideB);
dp[target] = modAdd(dp[target], modMul(structureWeight, ways[i]));
}
}
}
}
}
}
}
int answer = 0;
for (int used = 0; used <= n; ++used) {
answer = modAdd(answer, dp[index(used, 0, 0, 0, strideLen, strideA, strideB)]);
}
return answer;
}
private static int index(
int used,
int leftActive,
int rightActive,
int delta,
int strideLen,
int strideA,
int strideB) {
return used * strideLen + leftActive * strideA + rightActive * strideB + delta + DELTA_LIMIT;
}
}
private static void usage() {
System.err.println(
"Usage:\n"
+ " java Euler990 [n] [--skip-checkpoints] [--single-thread] [--threads=N]");
}
private static Options parseOptions(String[] args) {
Options options = new Options();
for (String arg : args) {
if ("--skip-checkpoints".equals(arg)) {
options.runCheckpoints = false;
continue;
}
if ("--single-thread".equals(arg)) {
options.allowMultithreading = false;
options.requestedThreads = 1;
continue;
}
if (arg.startsWith("--threads=")) {
options.requestedThreads = Integer.parseInt(arg.substring(10));
continue;
}
if (!arg.isEmpty() && arg.charAt(0) == '-') {
usage();
System.exit(1);
}
options.targetN = Integer.parseInt(arg);
}
if (options.targetN < 0 || options.targetN > TARGET_N) {
System.err.println("n must satisfy 0 <= n <= " + TARGET_N + ".");
System.exit(1);
}
return options;
}
private static void runCheckpoints(Solver solver) {
assert solver.solve(0) == 0;
assert solver.solve(1) == 0;
assert solver.solve(2) == 0;
assert solver.solve(3) == 9;
assert solver.solve(5) == 171;
assert solver.solve(7) == 4878;
}
public static void main(String[] args) {
Options options = parseOptions(args);
Solver solver = new Solver(options);
if (options.runCheckpoints) {
runCheckpoints(solver);
}
System.out.println(solver.solve(options.targetN));
}
}