Problem 294: Sum of Digits - Experience #23

Project 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

For a positive integer $k$, $d(k)$ is the sum of its decimal digits, for example $d(42)=4+2=6$. For a positive integer $n$, $S(n)$ is the number of positive integers $k\lt10^n$ that are divisible by $23$ and have $d(k)=23$. Given $S(9)=263626$ and $S(42)=6377168878570056$, we must find $S(11^{12})$ modulo $10^9$, that is, its remainder after division by $10^9$. ^[1]

The whole solution in one sentence: the remainders of the powers of $10$ modulo $23$ repeat every $22$ positions, so we sort the digit positions into $22$ groups in which only the sum of the digits matters, count the ways to reach each digit sum in a group by repeatedly squaring a small polynomial, and combine the groups in a table of digit sums and remainders.

Mathematical Approach

We first count digit by digit and then use the repeating place values to replace the $n$ positions by $22$ classes.

1. Digit strings and small cases

Write every $k\lt10^n$ with exactly $n$ digits, using leading zeros, with $d_0$ the units digit:

$$k=\sum_{i=0}^{n-1}d_i\,10^i,\qquad 0\le d_i\le9.$$

Leading zeros do not change the digit sum, and the string $00\dots0$, which is not a positive integer, has digit sum $0$. So $S(n)$ is the number of strings of $n$ digits with

$$\sum_{i=0}^{n-1}d_i=23,$$

$$\sum_{i=0}^{n-1}d_i\,10^i\equiv0\pmod{23},$$

where $a\equiv b\pmod{23}$ means that $a-b$ is a multiple of $23$. One or two digits add up to at most $18$, so $S(1)=S(2)=0$. Three digits with sum $23$ are an arrangement of $5,9,9$, of $6,8,9$, of $7,7,9$ or of $7,8,8$, which gives $3+6+3+3=15$ strings; none of them is a multiple of $23$, so $S(3)=0$. The smallest number with both properties is $1679=23\cdot73$, and a brute force over the multiples of $23$ gives $S(4)=20$, $S(5)=258$, $S(6)=2073$, $S(7)=12472$, $S(8)=61697$ and $S(9)=263626$, the value of the statement.

2. Counting digit by digit

Read the digits from $d_0$ upwards. After $d_0,\dots,d_{i-1}$ only two numbers matter: their sum $s$, where sums above $23$ can be dropped because a digit sum never decreases, and the remainder $t$ of $d_0+10\,d_1+\dots+10^{i-1}d_{i-1}$ modulo $23$. The digit $d$ at position $i$ changes this state by

$$\begin{aligned}s&\mapsto s+d,\\t&\mapsto(t+d\cdot10^i)\bmod23.\end{aligned}$$

If $C_i(s,t)$ is the number of strings of $i$ digits in the state $(s,t)$, then $C_0(0,0)=1$, each $C_{i+1}$ follows from $C_i$ by this rule, and $S(n)=C_n(23,0)$. This is dynamic programming: counts for longer strings are built from stored counts for shorter ones, and no string is ever listed. ^[2] With $24\cdot23=552$ states and at most $10$ digits each, a position costs a few thousand additions, and $S(9)$ and $S(42)$ appear at once. For $n=11^{12}=3138428376721$, however, the method would need more than three trillion steps.

3. The place values repeat every 22 positions

Since $23$ is a prime that does not divide $10$, Fermat's little theorem gives $10^{22}\equiv1\pmod{23}$. ^[3] Hence $10^{i+22}\equiv10^i$, and the place value of position $i$ modulo $23$ depends only on $c=i\bmod22$. The weights $w_c=10^c\bmod23$, found by multiplying by $10$ and reducing, are for $c=0,1,\dots,21$

$1,10,8,11,18,19,6,14,2,20,16,$
$22,13,15,12,5,4,17,9,21,3,7.$

Every number from $1$ to $22$ appears exactly once, so $22$ is the shortest period: $10$ is a primitive root modulo $23$, $23$ is a full reptend prime, and $1/23$ repeats the $22$ decimal digits $0434782608695652173913$, which form the number $(10^{22}-1)/23$. ^[4] The method uses only $10^{22}\equiv1$.

4. The idea in miniature: divisibility by 11

Modulo $11$ the place values repeat every $2$ positions, because $10\equiv-1\pmod{11}$: they are $1,-1,1,-1,\dots$, which explains the familiar rule that a number is divisible by $11$ exactly when its alternating digit sum is. ^[5] Let us count the multiples of $11$ below $10^4$ with digit sum $4$. Let $a$ be the sum of the digits in the even positions $0$ and $2$, and $b$ the sum in the odd positions $1$ and $3$. The number is congruent to $a-b$ modulo $11$, so we need $a+b=4$ and $a-b\equiv0\pmod{11}$. Since $a-b$ lies between $-4$ and $4$, this forces $a=b=2$. Two digits with sum $2$ are $02$, $11$ or $20$, so there are $3\cdot3=9$ such numbers:

$22$, $121$, $220$, $1012$, $1111$, $1210$, $2002$, $2101$ and $2200$.

We never looked at the four digits one by one, only at the sum of each group and the number of ways to reach it. For $23$ the same idea works with $22$ groups.

5. Classes of positions

Write $n=22q+r$ with $0\le r\le21$ and put position $i$ into the class $c=i\bmod22$. Class $c$ then has

$$u_c=\begin{cases}q+1,&c\lt r,\\q,&c\ge r\end{cases}$$

positions. For example, $n=25$ gives $q=1$ and $r=3$: the classes $0$, $1$ and $2$ hold the positions $\{0,22\}$, $\{1,23\}$ and $\{2,24\}$, the other classes one position each. For $n=11^{12}$ we get $q=142655835305$ and $r=11$, so the classes $0$ to $10$ have $q+1$ positions and the classes $11$ to $21$ have $q$. Let $k_c$ be the sum of the digits in class $c$, and $N_u(k)$ the number of strings of $u$ digits with digit sum $k$. An empty class, as for $n\lt22$, has $N_0(0)=1$ and $N_0(k)=0$ for $k\ge1$.

Claim. For every $22$-tuple of non-negative class sums $k_0,\dots,k_{21}$, multiply the $22$ counts $N_{u_c}(k_c)$. Then $S(n)$ is the sum of these products over the tuples with

$$\begin{aligned}\sum_{c=0}^{21}k_c&=23,\\\sum_{c=0}^{21}w_c\,k_c&\equiv0\pmod{23}.\end{aligned}$$

Proof. Every position of class $c$ has the place value $10^i\equiv w_c$, so

$$\sum_{i=0}^{n-1}d_i\,10^i\equiv\sum_{c=0}^{21}w_c\,k_c\pmod{23},$$

and the digit sum is $k_0+\dots+k_{21}$: both conditions of Section 1 depend only on the class sums. The classes split the positions, so a string of $n$ digits is the same as a choice of one digit string for each class. For fixed class sums these choices are independent, with $N_{u_c}(k_c)$ possibilities for class $c$, which gives the product, and every string is counted once, for its own class sums.

6. Counting digit strings with a polynomial

The numbers $N_u(k)$ are coefficients of a power:

$$N_u(k)=[x^k]\,\bigl(1+x+\dots+x^9\bigr)^u,$$

where $[x^k]$ means the coefficient of $x^k$. Multiplying out the $u$ factors means choosing a term $x^d$ from each factor, that is, a digit for each position, and the exponents add up to the digit sum. This is the basic idea of generating functions: multiplying the polynomials of independent parts combines their counts. ^[6] For $u=2$ the coefficients are $1,2,\dots,9,10,9,\dots,2,1$, so $N_2(k)=k+1$ for $k\le9$ and $N_2(k)=19-k$ for $10\le k\le18$; for instance $N_2(2)=3$ counts the pairs $02$, $11$, $20$ of Section 4.

Only the coefficients of $x^0,\dots,x^{23}$ are needed, because the class sums add up to $23$. Exponents never decrease in a product, so a term of degree above $23$ never contributes to a coefficient of degree at most $23$, and we may cut every polynomial off after $x^{23}$, even after each multiplication. A cut-off polynomial is a list of $24$ numbers, however large $u$ is.

For $u\ge1$ there is also a closed formula, the classical count for dice. ^[7] Without the bound $9$, the sequences of $u$ non-negative integers with sum $m$ correspond to the rows of $m$ marks and $u-1$ separators, so there are $M_u(m)=\binom{m+u-1}{u-1}$ of them. Inclusion–exclusion over the entries that are at least $10$, subtracting $10$ from each of $j$ chosen entries, then gives

$$N_u(k)=\sum_{j=0}^{J}(-1)^j\binom{u}{j}M_u(k-10j)$$

where $J=\lfloor k/10\rfloor$, so $J\le2$ for $k\le23$. For example, two digits with sum $10$ give $M_2(10)-2\,M_2(0)$${}=11-2=9$. The programs do not use the formula. Its binomial coefficients are ratios of factorials, and remainders modulo $10^9$ cannot be divided by factorials, which contain the factors $2$ and $5$; the simplest way is to evaluate the formula exactly and reduce afterwards, as in the check below, while the polynomial powers need only additions and multiplications of remainders.

7. Huge powers by repeated squaring

The classes have only the sizes $q$ and $q+1$, so two cut-off powers of $B=1+x+\dots+x^9$ are needed, with $q$ near $1.4\cdot10^{11}$. Exponentiation by squaring reads the binary digits of the exponent from the lowest one up, squares $B$ to $B^2,B^4,B^8,\dots$ and multiplies the result by the powers that belong to the binary digits $1$. ^[8] For example $9=1001_2$, so $B^9=B\cdot B^8$, with three squarings and two multiplications into the result. The exponents $q=142655835305$ and $q+1$ both have $38$ binary digits, $19$ of them $1$, so each power costs $37$ squarings and $19$ further multiplications. Multiplying two cut-off polynomials $\sum a_ix^i$ and $\sum b_jx^j$ adds $a_ib_j$ to the coefficient of $x^{i+j}$ for the $300$ pairs with $i+j\le23$; the programs skip zero coefficients, so at most $300$ products are formed.

8. Combining the classes

The sum of Section 5 is built up class by class, like the count of Section 2 with a class in the role of a digit. Let $D_c(s,t)$ be the number of ways to fill the first $c$ classes with digits so that their sum is $s$ and their contribution to the number has the remainder $t$ modulo $23$; by the proof of Section 5, it is the sum of the products $N_{u_0}(k_0)\cdots N_{u_{c-1}}(k_{c-1})$ over the matching class sums. We start with $D_0(0,0)=1$ and all other entries $0$. To add class $c$, we start a fresh table of zeros, and every state $(s,t)$ and every class sum $k$ with $s+k\le23$ contribute $D_c(s,t)\,N_{u_c}(k)$ to the state

$$\bigl(s+k,\;(t+w_c\,k)\bmod23\bigr).$$

After the $22$ classes, $S(n)=D_{22}(23,0)$ by the claim of Section 5.

The example of Section 4 in this form has two classes with the weights $1$ and $10\equiv-1\pmod{11}$, two positions each, and the cut-off polynomial $1+2x+3x^2+4x^3+5x^4$. The first class puts $N_2(k)$ into the state $(k,k)$, for instance $5$ into $(4,4)$. After the second class the states with digit sum $4$ hold

remainder $0$: $3\cdot3=9$, from $a=b=2$;
remainder $2$: $4\cdot2=8$, from $a=3$ and $b=1$;
remainder $4$: $5\cdot1=5$, from $a=4$ and $b=0$;
remainder $7\equiv-4$: $1\cdot5=5$, from $a=0$ and $b=4$;
remainder $9\equiv-2$: $2\cdot4=8$, from $a=1$ and $b=3$.

The entry for remainder $0$ is the count of Section 4, and the five entries add up to $35$, the number of all strings of four digits with digit sum $4$, as a brute force confirms.

9. Only the last nine digits

The exact value of $S(11^{12})$ has $264$ digits, as an exact computation with big integers shows, but only its remainder modulo $10^9$ is asked for. Every step above is an addition or a multiplication, and if $a\equiv a'$ and $b\equiv b'$, then $a+b\equiv a'+b'$ and $ab\equiv a'b'$, because $ab-a'b'=a(b-b')+(a-a')b'$. So every number may be replaced by its remainder modulo $10^9$ as soon as it appears; no division occurs anywhere. A product of two remainders is below $10^{18}\lt2^{63}$, so it fits into a signed 64-bit integer.

How the Code Works

The C++ program ^[9] starts in main. parse_arguments accepts --skip-checkpoints, --n= and --mod=, with the defaults kDefaultN $=11^{12}$ and kDefaultMod $=10^9$; any other argument, or the modulus $0$, gives exit code $1$. Unless skipped, run_checkpoints runs, and a failure gives exit code $2$. Then main prints count_via_classes_mod(n, mod).

A Coeff is an std::array of $24$ unsigned 64-bit numbers, a cut-off polynomial of Section 6. poly_mul_mod multiplies two of them over all $i+j\le23$, forming each product as an unsigned __int128 ^[10] and reducing it modulo mod. poly_pow_digit_sum runs the method of Section 7, starting from base $=1+x+\dots+x^9$ and result $=1$, and skips the squaring after the last bit.

count_via_classes_mod computes $q$, $r$, the powers c_q and c_q1 and the array weights of the $w_c$. Its tables dp and next have $24\times23$ entries, indexed by digit sum and remainder, and dp[0][0] starts at $1$. Class cls uses c_q1 if cls < r and c_q otherwise; the rule of Section 8 fills the cleared next from the non-zero entries of dp, again with 128-bit products, and next then replaces dp. The result is dp[23][0]. Every reduced number is below mod, and the sum of two of them is formed in 64 bits, so --mod= is safe for every modulus up to $2^{63}$. Larger moduli are accepted but can overflow these sums, and parse_u64_after_prefix lets decimal values above $2^{64}-1$ wrap around.

The checkpoints use Section 2: direct_count_exact keeps exact counts in unsigned __int128, and direct_count_mod reduces them modulo a given number. run_checkpoints requires direct_count_exact(9) $=263626$ and direct_count_exact(42) $=6377168878570056$, the values of the statement found without classes, and requires direct_count_mod and count_via_classes_mod to agree for $n=200$ modulo $10^9$, a case with $q=9$ and $r=2$.

The Python and Java files contain only the fast method, with the same three functions, called polyMulMod, polyPowDigitSum and countViaClassesMod in Java, and a solve that uses $n=11^{12}$ and the modulus $10^9$, without options or checkpoints. Python's integers are unbounded. Java multiplies in long and switches to BigInteger when a product comes out negative; with the modulus $10^9$ every product of two remainders is below $10^{18}\lt2^{63}$, so this branch never runs, and for larger moduli the sign test would not catch every overflow.

Complexity and Verification

The two powers take $2\cdot(37+19)=112$ multiplications with at most $300$ products each, so this part grows like $T^2\log n$ for the digit sum $T=23$. The class table needs at most $22\cdot23\cdot300=151800$ updates, for $22$ classes, $23$ remainders and $300$ pairs $(s,k)$, whatever $n$ is. The memory does not depend on $n$: two tables of $552$ numbers, the $22$ weights and a few lists of $24$ coefficients. The checkpoints cost more: the count of Section 2 makes at most $23\cdot195=4485$ updates per position, for the $195$ pairs $(s,d)$ with $s+d\le23$ and $d\le9$, over $9+42+200=251$ positions, and the class method runs once more for $n=200$. On an Apple M4 Max the C++ program runs in about $4$ milliseconds with its checkpoints, the Python and Java programs in about $30$, mostly start-up time.

Besides the checkpoints, the following checks were run for this article; they are not part of the published solutions.

A brute force over all multiples of $23$ below $10^9$ gives the values $S(1)$ to $S(9)$ of Section 1.
The class method with exact big integers and the count of Section 2 agree for every $n\le300$, including $S(42)$, and the Python count_via_classes_mod agrees with the exact counts modulo $10^9$ for every $n\le400$.
For the exponents $q$ and $q+1$ of $n=11^{12}$, the coefficients from poly_pow_digit_sum agree with the formula of Section 6, evaluated exactly and reduced modulo $10^9$.
A decomposition without classes cuts the positions into $q$ blocks of $22$ consecutive positions and $r$ left over, describes a block by a polynomial in $x$ for the digit sum and $y$ for the remainder, with $y^{23}=1$, and raises it to the power $q$. It agrees with the exact counts for $n\le120$ and with the C++ result for $n=11^{12}$.
The C++ program, built with clang++ and with g++, and the Python and Java programs print the same value, and so does the exact computation of Section 9.

Footnotes and References

The references supply the background; the reduction to class sums, the cut-off powers, the class table and all counts are derived or computed in this article.

Project Euler 294 — Sum of Digits - Experience #23. The official statement defines d(k) as the decimal digit sum and S(n) as the number of positive integers k below 10ⁿ divisible by 23 with d(k) = 23, gives S(9) = 263626 and S(42) = 6377168878570056, and asks for S(11¹²) modulo 10⁹.
Wikipedia — Dynamic programming. A problem is broken down into simpler sub-problems whose solutions are stored and reused, so that each sub-problem is solved once.
Wikipedia — Fermat's little theorem. If p is a prime and a is not divisible by p, then a^(p−1) ≡ 1 (mod p).
Wikipedia — Full reptend prime. In base 10 these are the primes p for which 10 is a primitive root modulo p, so that 1/p has a decimal period of p − 1 digits; 23 is in the list, OEIS A001913.
Wikipedia — Divisibility rule. A number is divisible by 11 exactly when the alternating sum of its digits is, because 10 ≡ −1 (mod 11).
Wikipedia — Generating function. A sequence is stored as the coefficients of a power series, and multiplying two ordinary generating functions gives the convolution of their sequences.
MathWorld — Dice. The number of ways to roll the total p with n dice of s faces is the coefficient of x^p in (x + x² + … + x^s)^n, the sum over 0 ≤ i ≤ ⌊(p − n)/s⌋ of (−1)^i C(n, i) C(p − si − 1, n − 1). For s = 10, n = u and p = k + u it is the formula of Section 6.
Wikipedia — Exponentiation by squaring. A power xⁿ needs about log₂ n squarings and at most as many further multiplications, in any semigroup, for example of polynomials.
C++ — Euler294.cpp. The linked immutable revision is the implementation record for this article.
GCC manual — 128-bit Integers. GCC supports __int128 and unsigned __int128 on targets with an integer mode wide enough to hold 128 bits.

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