Problem 294: Sum of Digits - Experience #23
View on Project EulerProject Euler Problem 294 Solution
EulerSolve provides an optimized solution for Project Euler Problem 294, Sum of Digits - Experience #23, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Let \(S(n)\) denote the number of integers $$0 \le x < 10^n$$ whose decimal digit sum is \(23\) and which are divisible by \(23\). Leading zeros are therefore allowed in the \(n\)-digit representation, which matches the code. The Project Euler target is \(S(11^{12}) \bmod 10^9\). Mathematical Approach 1. Writing the two conditions in digit form Write $$x = \sum_{i=0}^{n-1} d_i 10^i,\qquad 0 \le d_i \le 9,$$ where \(d_0\) is the units digit. Then the problem imposes two simultaneous constraints: $$\sum_{i=0}^{n-1} d_i = 23,$$ $$\sum_{i=0}^{n-1} d_i 10^i \equiv 0 \pmod{23}.$$ A direct dynamic program over all \(n\) positions would be impossible for \(n = 11^{12}\), so the key is to compress positions that behave identically modulo \(23\). 2. Why positions fall into 22 residue classes Modulo \(23\), powers of \(10\) repeat with period \(22\), because $$10^{22} \equiv 1 \pmod{23}.$$ Therefore, if two positions satisfy \(i \equiv j \pmod{22}\), then $$10^i \equiv 10^j \pmod{23}.$$ So all digits sitting in the same residue class modulo \(22\) contribute to divisibility by the same weight. Define $$w_c = 10^c \bmod 23,\qquad c=0,1,\ldots,21.$$ If \(n = 22q + r\) with \(0 \le r < 22\), then class \(c\) contains $$u_c = \begin{cases} q+1, & c < r, \\ q, & c \ge r....
Detailed mathematical approach
Problem Summary
Let \(S(n)\) denote the number of integers
$$0 \le x < 10^n$$
whose decimal digit sum is \(23\) and which are divisible by \(23\). Leading zeros are therefore allowed in the \(n\)-digit representation, which matches the code. The Project Euler target is \(S(11^{12}) \bmod 10^9\).
Mathematical Approach
1. Writing the two conditions in digit form
Write
$$x = \sum_{i=0}^{n-1} d_i 10^i,\qquad 0 \le d_i \le 9,$$
where \(d_0\) is the units digit. Then the problem imposes two simultaneous constraints:
$$\sum_{i=0}^{n-1} d_i = 23,$$
$$\sum_{i=0}^{n-1} d_i 10^i \equiv 0 \pmod{23}.$$
A direct dynamic program over all \(n\) positions would be impossible for \(n = 11^{12}\), so the key is to compress positions that behave identically modulo \(23\).
2. Why positions fall into 22 residue classes
Modulo \(23\), powers of \(10\) repeat with period \(22\), because
$$10^{22} \equiv 1 \pmod{23}.$$
Therefore, if two positions satisfy \(i \equiv j \pmod{22}\), then
$$10^i \equiv 10^j \pmod{23}.$$
So all digits sitting in the same residue class modulo \(22\) contribute to divisibility by the same weight. Define
$$w_c = 10^c \bmod 23,\qquad c=0,1,\ldots,21.$$
If \(n = 22q + r\) with \(0 \le r < 22\), then class \(c\) contains
$$u_c = \begin{cases} q+1, & c < r, \\ q, & c \ge r. \end{cases}$$
For example, if \(n=25\), then \(q=1\), \(r=3\), so classes \(0,1,2\) contain two positions each, while the other \(19\) classes contain one position each.
3. Reducing each class to its digit sum
Inside one class \(c\), suppose the digits in that class add up to \(k\). Then their total contribution modulo \(23\) is simply
$$k\,w_c \pmod{23},$$
because every position in that class carries the same weight \(w_c\). So we do not need the full digit pattern of the class; we only need to know how many assignments produce class sum \(k\).
For a class of size \(u\), the generating polynomial is
$$G_u(x) = (1+x+x^2+\cdots+x^9)^u.$$
The coefficient \([x^k]G_u(x)\) is exactly the number of ways to place \(u\) digits from \(0\) to \(9\) whose sum is \(k\).
Because the global digit sum is only \(23\), coefficients above degree \(23\) can never matter. So the code truncates every polynomial to degree \(23\), turning the huge exponent \(u\) into a manageable binary-exponentiation problem.
4. Why only two class polynomials are needed
Every class size is either \(q\) or \(q+1\). Therefore we only need the two truncated polynomials
$$G_q(x) = (1+x+\cdots+x^9)^q,$$
$$G_{q+1}(x) = (1+x+\cdots+x^9)^{q+1},$$
both computed modulo \(10^9\). The function poly_pow_digit_sum does this by repeated squaring, and the multiplication routine discards all terms above \(x^{23}\).
5. DP over classes: total digit sum and remainder mod 23
After compressing each residue class into “choose a class sum \(k\) with multiplicity \([x^k]G_{u_c}(x)\)”, the problem becomes a 22-step dynamic program. Let
$$DP[s][t]$$
be the number of ways, after processing some classes, to obtain total digit sum \(s\) and remainder \(t \pmod{23}\).
When processing class \(c\), choosing class sum \(k\) produces the transition
$$s' = s+k,\qquad t' = (t + w_c k) \bmod 23.$$
The multiplicity of this transition is \([x^k]G_{u_c}(x)\). After all \(22\) classes have been processed, the desired count is
$$DP[23][0].$$
This is the whole reason the huge value of \(n\) becomes harmless: the algorithm never iterates over all positions, only over the \(22\) residue classes.
6. Checks and interpretation
The implementation verifies itself with three checkpoints:
$$S(9)=263626,$$
$$S(42)=6377168878570056,$$
and a consistency check that the compressed class DP matches a direct position-by-position DP for
$$n=200.$$
The direct DP is only feasible for moderate \(n\), but it is an excellent validation tool: it confirms that grouping positions by residue modulo \(22\) loses no information.
How the Code Works
poly_mul_mod multiplies two digit-sum polynomials and truncates to degree \(23\). poly_pow_digit_sum raises \(1+x+\cdots+x^9\) to exponent \(q\) or \(q+1\) by binary exponentiation. count_via_classes_mod builds the \(22\) modular weights \(w_c\), runs the DP over \((\text{digit sum}, \text{remainder})\), and returns \(DP[23][0]\). The functions direct_count_exact and direct_count_mod are only for checkpoints on smaller \(n\).
Complexity Analysis
Polynomial multiplication is on vectors of length \(24\), so one multiplication costs \(O(23^2)\). Binary exponentiation therefore costs
$$O(23^2 \log n).$$
The class DP visits \(22\) classes, \(24\) digit-sum states, \(23\) remainders, and up to \(24\) class-sum choices, so its cost is
$$O(22 \cdot 23^3),$$
which is effectively constant for this fixed target sum. Memory usage is \(O(23^2)\).
Further Reading
- Problem page: https://projecteuler.net/problem=294
- Generating functions in combinatorics: https://en.wikipedia.org/wiki/Generating_function
- Modular arithmetic and multiplicative order: https://en.wikipedia.org/wiki/Multiplicative_order
Problem 294 source code
C++
#include <array>
#include <cstdint>
#include <iostream>
#include <string>
#include <algorithm>
#include <functional>
namespace {
using u64 = std::uint64_t;
using u128 = unsigned __int128;
constexpr int kTargetSum = 23;
constexpr int kMod23 = 23;
constexpr u64 kDefaultMod = 1000000000ULL;
constexpr u64 kDefaultN = 3138428376721ULL; // 11^12
struct Options {
u64 n = kDefaultN;
u64 mod = kDefaultMod;
bool run_checkpoints = true;
};
bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& value) {
if (arg.rfind(prefix, 0U) != 0U) {
return false;
}
const std::string tail = arg.substr(prefix.size());
if (tail.empty()) {
return false;
}
u64 parsed = 0;
for (char c : tail) {
if (c < '0' || c > '9') {
return false;
}
parsed = parsed * 10 + static_cast<u64>(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_u64_after_prefix(arg, "--n=", options.n)) {
continue;
}
if (parse_u64_after_prefix(arg, "--mod=", options.mod)) {
continue;
}
std::cerr << "Unknown argument: " << arg << '\n';
return false;
}
return options.mod > 0;
}
std::string to_string_u128(u128 value) {
if (value == 0) {
return "0";
}
std::string out;
while (value > 0) {
const int digit = static_cast<int>(value % 10);
out.push_back(static_cast<char>('0' + digit));
value /= 10;
}
std::reverse(out.begin(), out.end());
return out;
}
using Coeff = std::array<u64, kTargetSum + 1>;
Coeff poly_mul_mod(const Coeff& a, const Coeff& b, const u64 mod) {
Coeff c{};
c.fill(0);
for (int i = 0; i <= kTargetSum; ++i) {
if (a[static_cast<std::size_t>(i)] == 0) {
continue;
}
for (int j = 0; i + j <= kTargetSum; ++j) {
if (b[static_cast<std::size_t>(j)] == 0) {
continue;
}
const u128 add = static_cast<u128>(a[static_cast<std::size_t>(i)]) *
static_cast<u128>(b[static_cast<std::size_t>(j)]);
c[static_cast<std::size_t>(i + j)] =
(c[static_cast<std::size_t>(i + j)] + static_cast<u64>(add % mod)) % mod;
}
}
return c;
}
Coeff poly_pow_digit_sum(const u64 exponent, const u64 mod) {
Coeff base{};
base.fill(0);
for (int d = 0; d <= 9; ++d) {
base[static_cast<std::size_t>(d)] = 1;
}
Coeff result{};
result.fill(0);
result[0] = 1;
u64 e = exponent;
while (e > 0) {
if ((e & 1ULL) != 0ULL) {
result = poly_mul_mod(result, base, mod);
}
e >>= 1U;
if (e > 0) {
base = poly_mul_mod(base, base, mod);
}
}
return result;
}
u64 count_via_classes_mod(const u64 n, const u64 mod) {
const u64 q = n / 22ULL;
const int r = static_cast<int>(n % 22ULL);
const Coeff c_q = poly_pow_digit_sum(q, mod);
const Coeff c_q1 = poly_pow_digit_sum(q + 1ULL, mod);
std::array<int, 22> weights{};
int cur = 1;
for (int i = 0; i < 22; ++i) {
weights[static_cast<std::size_t>(i)] = cur;
cur = (cur * 10) % kMod23;
}
std::array<std::array<u64, kMod23>, kTargetSum + 1> dp{};
std::array<std::array<u64, kMod23>, kTargetSum + 1> next{};
for (auto& row : dp) {
row.fill(0);
}
dp[0][0] = 1;
for (int cls = 0; cls < 22; ++cls) {
const Coeff& coeff = (cls < r) ? c_q1 : c_q;
for (auto& row : next) {
row.fill(0);
}
for (int s = 0; s <= kTargetSum; ++s) {
for (int rem = 0; rem < kMod23; ++rem) {
const u64 cur_count = dp[static_cast<std::size_t>(s)][static_cast<std::size_t>(rem)];
if (cur_count == 0) {
continue;
}
for (int t = 0; s + t <= kTargetSum; ++t) {
const u64 ways = coeff[static_cast<std::size_t>(t)];
if (ways == 0) {
continue;
}
const int rem2 = (rem + weights[static_cast<std::size_t>(cls)] * t) % kMod23;
const u128 add = static_cast<u128>(cur_count) * static_cast<u128>(ways);
next[static_cast<std::size_t>(s + t)][static_cast<std::size_t>(rem2)] =
(next[static_cast<std::size_t>(s + t)][static_cast<std::size_t>(rem2)] +
static_cast<u64>(add % mod)) % mod;
}
}
}
dp = next;
}
return dp[kTargetSum][0];
}
u128 direct_count_exact(const int n) {
std::array<std::array<u128, kMod23>, kTargetSum + 1> dp{};
std::array<std::array<u128, kMod23>, kTargetSum + 1> next{};
for (auto& row : dp) {
row.fill(0);
}
dp[0][0] = 1;
int pow10 = 1;
for (int pos = 0; pos < n; ++pos) {
for (auto& row : next) {
row.fill(0);
}
for (int s = 0; s <= kTargetSum; ++s) {
for (int rem = 0; rem < kMod23; ++rem) {
const u128 cur_count = dp[static_cast<std::size_t>(s)][static_cast<std::size_t>(rem)];
if (cur_count == 0) {
continue;
}
for (int d = 0; s + d <= kTargetSum && d <= 9; ++d) {
const int rem2 = (rem + d * pow10) % kMod23;
next[static_cast<std::size_t>(s + d)][static_cast<std::size_t>(rem2)] += cur_count;
}
}
}
dp = next;
pow10 = (pow10 * 10) % kMod23;
}
return dp[kTargetSum][0];
}
u64 direct_count_mod(const int n, const u64 mod) {
std::array<std::array<u64, kMod23>, kTargetSum + 1> dp{};
std::array<std::array<u64, kMod23>, kTargetSum + 1> next{};
for (auto& row : dp) {
row.fill(0);
}
dp[0][0] = 1;
int pow10 = 1;
for (int pos = 0; pos < n; ++pos) {
for (auto& row : next) {
row.fill(0);
}
for (int s = 0; s <= kTargetSum; ++s) {
for (int rem = 0; rem < kMod23; ++rem) {
const u64 cur_count = dp[static_cast<std::size_t>(s)][static_cast<std::size_t>(rem)];
if (cur_count == 0) {
continue;
}
for (int d = 0; s + d <= kTargetSum && d <= 9; ++d) {
const int rem2 = (rem + d * pow10) % kMod23;
next[static_cast<std::size_t>(s + d)][static_cast<std::size_t>(rem2)] += cur_count;
if (next[static_cast<std::size_t>(s + d)][static_cast<std::size_t>(rem2)] >= mod) {
next[static_cast<std::size_t>(s + d)][static_cast<std::size_t>(rem2)] %= mod;
}
}
}
}
dp = next;
pow10 = (pow10 * 10) % kMod23;
}
return dp[kTargetSum][0] % mod;
}
bool run_checkpoints() {
const u128 s9 = direct_count_exact(9);
if (s9 != static_cast<u128>(263626ULL)) {
std::cerr << "Checkpoint failed for S(9): got " << to_string_u128(s9) << '\n';
return false;
}
const u128 s42 = direct_count_exact(42);
if (s42 != static_cast<u128>(6377168878570056ULL)) {
std::cerr << "Checkpoint failed for S(42): got " << to_string_u128(s42) << '\n';
return false;
}
const u64 direct200 = direct_count_mod(200, kDefaultMod);
const u64 class200 = count_via_classes_mod(200, kDefaultMod);
if (direct200 != class200) {
std::cerr << "Class-DP mismatch at n=200: direct=" << direct200
<< ", class=" << class200 << '\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;
}
const u64 answer = count_via_classes_mod(options.n, options.mod);
std::cout << answer << '\n';
return 0;
}
Python
def poly_mul_mod(a, b, mod):
c = [0] * 24
for i in range(24):
if a[i] == 0:
continue
for j in range(24 - i):
if b[j] == 0:
continue
c[i + j] = (c[i + j] + a[i] * b[j]) % mod
return c
def poly_pow_digit_sum(exponent, mod):
base = [0] * 24
for d in range(10):
base[d] = 1
result = [0] * 24
result[0] = 1
e = exponent
while e > 0:
if e & 1:
result = poly_mul_mod(result, base, mod)
e >>= 1
if e > 0:
base = poly_mul_mod(base, base, mod)
return result
def count_via_classes_mod(n, mod):
q = n // 22
r = int(n % 22)
c_q = poly_pow_digit_sum(q, mod)
c_q1 = poly_pow_digit_sum(q + 1, mod)
weights = [0] * 22
cur = 1
for i in range(22):
weights[i] = cur
cur = (cur * 10) % 23
dp = [[0] * 23 for _ in range(24)]
dp[0][0] = 1
for cls in range(22):
coeff = c_q1 if cls < r else c_q
nxt = [[0] * 23 for _ in range(24)]
for s in range(24):
for rem in range(23):
cur_count = dp[s][rem]
if cur_count == 0:
continue
w = weights[cls]
for t in range(24 - s):
ways = coeff[t]
if ways == 0:
continue
rem2 = (rem + w * t) % 23
nxt[s + t][rem2] = (nxt[s + t][rem2] + cur_count * ways) % mod
dp = nxt
return dp[23][0]
def solve(n=3138428376721, mod=1000000000):
return str(count_via_classes_mod(n, mod))
if __name__ == '__main__':
print(solve())
Java
public class Euler294 {
static final int TARGET_SUM = 23;
static final int MOD_23 = 23;
static long[] polyMulMod(long[] a, long[] b, long mod) {
long[] c = new long[TARGET_SUM + 1];
for (int i = 0; i <= TARGET_SUM; ++i) {
if (a[i] == 0)
continue;
for (int j = 0; i + j <= TARGET_SUM; ++j) {
if (b[j] == 0)
continue;
long term = (a[i] % mod) * (b[j] % mod);
if (term < 0) { // Overflow prevention fallback
java.math.BigInteger bTerm = java.math.BigInteger.valueOf(a[i])
.multiply(java.math.BigInteger.valueOf(b[j]))
.mod(java.math.BigInteger.valueOf(mod));
c[i + j] = (c[i + j] + bTerm.longValue()) % mod;
} else {
c[i + j] = (c[i + j] + term % mod) % mod;
}
}
}
return c;
}
static long[] polyPowDigitSum(long exponent, long mod) {
long[] base = new long[TARGET_SUM + 1];
for (int d = 0; d <= 9; ++d) {
base[d] = 1;
}
long[] result = new long[TARGET_SUM + 1];
result[0] = 1;
long e = exponent;
while (e > 0) {
if ((e & 1) != 0) {
result = polyMulMod(result, base, mod);
}
e >>= 1;
if (e > 0) {
base = polyMulMod(base, base, mod);
}
}
return result;
}
static long countViaClassesMod(long n, long mod) {
long q = n / 22L;
int r = (int) (n % 22L);
long[] cQ = polyPowDigitSum(q, mod);
long[] cQ1 = polyPowDigitSum(q + 1L, mod);
int[] weights = new int[22];
int cur = 1;
for (int i = 0; i < 22; ++i) {
weights[i] = cur;
cur = (cur * 10) % MOD_23;
}
long[][] dp = new long[TARGET_SUM + 1][MOD_23];
dp[0][0] = 1;
for (int cls = 0; cls < 22; ++cls) {
long[] coeff = (cls < r) ? cQ1 : cQ;
long[][] next = new long[TARGET_SUM + 1][MOD_23];
for (int s = 0; s <= TARGET_SUM; ++s) {
for (int rem = 0; rem < MOD_23; ++rem) {
long curCount = dp[s][rem];
if (curCount == 0)
continue;
for (int t = 0; s + t <= TARGET_SUM; ++t) {
long ways = coeff[t];
if (ways == 0)
continue;
int rem2 = (rem + weights[cls] * t) % MOD_23;
long term = (curCount % mod) * (ways % mod);
if (term < 0) {
java.math.BigInteger bTerm = java.math.BigInteger.valueOf(curCount)
.multiply(java.math.BigInteger.valueOf(ways))
.mod(java.math.BigInteger.valueOf(mod));
next[s + t][rem2] = (next[s + t][rem2] + bTerm.longValue()) % mod;
} else {
next[s + t][rem2] = (next[s + t][rem2] + term % mod) % mod;
}
}
}
}
dp = next;
}
return dp[TARGET_SUM][0];
}
public static String solve() {
return String.valueOf(countViaClassesMod(3138428376721L, 1000000000L));
}
public static void main(String[] args) {
System.out.println(solve());
}
}