Problem 377: Sum of Digits - Experience #13

Project Euler Problem 377 Solution

EulerSolve provides an optimized solution for Project Euler Problem 377, Sum of Digits - Experience #13, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For a positive integer $n$, let $F(n)=f(n)$ be the sum of all positive decimal integers that contain no zero digit and whose digit sum is exactly $n$. Let $A(n)$ denote how many such integers exist. The problem asks for $$\sum_{k=1}^{17} F(13^k)\pmod{10^9}.$$ A direct enumeration is impossible because $13^{17}$ is enormous, so the implementation turns the digit-building process into a fixed linear recurrence and then evaluates that recurrence with matrix exponentiation. Mathematical Approach Step 1: Count the admissible numbers Every valid number with digit sum $n$ ends in exactly one digit $d\in\{1,\dots,9\}$. After removing that last digit, the remaining prefix must still avoid zeros and must have digit sum $n-d$. Therefore $$A(n)=\sum_{d=1}^{9} A(n-d).$$ To make the recurrence uniform, we use the standard empty-prefix convention $$A(0)=1,\qquad A(n)=0 \text{ for } n\lt 0.$$ The value $A(0)=1$ does not represent a positive integer. It represents the single empty prefix that allows one-digit numbers to be generated correctly. Step 2: Sum the values of the numbers Now let $F(n)$ be the sum of all valid integers with digit sum $n$. Fix the last digit $d$. If a prefix value is $x$, appending $d$ on the right creates the new number $10x+d$. Summing over all prefixes of digit sum $n-d$ gives the contribution $10F(n-d)+d\,A(n-d)$....

Detailed mathematical approach

Problem Summary

Call a positive integer zeroless when none of its decimal digits is $0$, and let $f(n)$ be the sum of all zeroless integers whose digits add up to $n$. The statement lists the $16$ zeroless integers of digital sum $5$ — $5$, $14$, $23$, $32$, $41$, $113$, $122$, $131$, $212$, $221$, $311$, $1112$, $1121$, $1211$, $2111$ and $11111$ — and gives $f(5)=17891$. We must evaluate

$$S=\sum_{i=1}^{17}f(13^i)$$

and report the last nine digits of $S$. ^[1]

The whole solution in one sentence: peeling the last digit off a zeroless integer turns both the count of these integers and their sum into recurrences that look back nine steps, so one fixed $18\times18$ matrix advances both, and one shared table of its powers, built by repeated squaring modulo $10^9$, answers all seventeen queries.

Mathematical Approach

1. Zeroless integers are compositions

A zeroless integer is nothing but its digit string $d_1d_2\ldots d_k$ with every $d_j$ between $1$ and $9$, of digit sum $d_1+\cdots+d_k$. An ordered list of positive integers with a given total is a composition of that total, ^[2] so the zeroless integers of digit sum $n$ are exactly the compositions of $n$ into parts at most $9$ — one-to-one, since a leading digit is never $0$ and distinct strings therefore denote distinct integers.

This explains the $16$ of the statement: there are $2^{n-1}$ compositions of $n$ in all, one per choice of cutting or joining each of the $n-1$ gaps between $n$ ones, ^[2] and for $n\le9$ no part can exceed $9$ anyway, so the count at $n=5$ is $2^4=16$. The bound first bites at $n=10$, where the composition $(10)$ must be dropped, leaving $2^9-1=511$.

2. Counting them: a nine-step recurrence

Let $a_n$ be the number of digit strings with digits in $\{1,\ldots,9\}$ and digit sum $n$. The last digit is some $d$ between $1$ and $9$, and removing it leaves a shorter string of digit sum $n-d$, so

$$a_n=\sum_{d=1}^{9}a_{n-d},\qquad a_0=1,$$

with $a_n=0$ for $n\lt0$. The value $a_0=1$ counts the empty string: not a positive integer, but what lets a one-digit integer appear as an empty prefix followed by one digit. For $n\ge1$, $a_n$ is the honest count. The first values are

$$\begin{aligned}a_1,\ldots,a_4&=1,2,4,8,\\ a_5,\ldots,a_8&=16,32,64,128,\\ a_9,\ldots,a_{12}&=256,511,1021,2040,\end{aligned}$$

powers of two up to $a_9$, then slower growth once long parts are forbidden. This nine-step Fibonacci rule is tabulated with the zeroless-integer meaning in the OEIS. ^[3]

3. Summing them: peel off the last digit

Let $f_n$ be the sum of the values of those same strings, with $f_0=0$ and $f_n=0$ for $n\lt0$. Write a string with last digit $d$ as $10m+d$, $m$ being the value of the prefix. Summing $10m+d$ over the $a_{n-d}$ prefixes of digit sum $n-d$ gives $10f_{n-d}+d\,a_{n-d}$, and adding the nine choices of $d$ yields

$$f_n=\sum_{d=1}^{9}\bigl(10f_{n-d}+d\,a_{n-d}\bigr).$$

The coupling runs one way only: $f$ needs $a$, but $a$ never needs $f$. Every nonempty string has exactly one last digit and one prefix, so the nine groups are disjoint and exhaust everything, and $f_n=f(n)$ for $n\ge1$. For $1\le n\le9$ the term $d=n$ contributes $10f_0+n\,a_0=n$, which is the one-digit integer $n$; for $n\ge10$ no one-digit integer occurs.

4. A worked example: $f(5)$ by hand

Start from $a_0,\ldots,a_4=1,1,2,4,8$ and $f_0,\ldots,f_4=0,1,13,147,1625$. The five terms of the recurrence for $n=5$ are the sums of the sixteen integers grouped by their last digit:

$d=1$: $10\cdot1625+1\cdot8=16258$;
$d=2$: $10\cdot147+2\cdot4=1478$;
$d=3$: $10\cdot13+3\cdot2=136$;
$d=4$: $10\cdot1+4\cdot1=14$;
$d=5$: $10\cdot0+5\cdot1=5$.

Each value can be read off the list of sixteen; the eight integers ending in $1$ do add up to $16258$,

$$\begin{aligned}&41+131+221+311\\&{}+1121+1211+2111+11111,\end{aligned}$$

and likewise $32+122+212+1112=1478$ and $23+113=136$. The total $16258+1478+136+14+5=17891$ is the checkpoint of the statement, and the same recurrence gives $f_6=196833$, $f_9=261994131$ and $f_{13}=3835856884757$.

5. A generating function, and the order of the recurrence

Where do these recurrences come from? Put

$$P(x)=\sum_{d=1}^{9}x^d,\qquad D(x)=\sum_{d=1}^{9}d\,x^d.$$

Multiplying polynomials adds exponents, so the coefficient of $x^n$ in $P(x)^k$ counts the digit strings of length $k$ and digit sum $n$: in $P(x)^2$ it is $2$ at $x^3$, for $12$ and $21$. Writing $A(x)=\sum_n a_nx^n$ and $F(x)=\sum_n f_nx^n$, summing $P(x)^k$ over all $k\ge0$ gives $A(x)=1/\bigl(1-P(x)\bigr)$. ^[2] For $F$, sum over pairs (integer, digit position): a digit $d$ with $t$ digits to its right contributes $d\cdot10^{t}$. Its prefix is an arbitrary string, the digit carries weight $d$, and a suffix of $t$ parts carries $10^{t}$, so the factors are $1/(1-P)$, $D$ and $\sum_{t\ge0}(10P)^t$:

$$F(x)=\frac{D(x)}{\bigl(1-P(x)\bigr)\bigl(1-10P(x)\bigr)}.$$

Expanded as power series, $A$ and $F$ reproduce $a_n$ and $f_n$ up to $n=40$, which we checked term by term. Put $Q=1-11P+10P^2$, of degree $18$, and equate coefficients of $x^n$ in $QF=D$. ^[6] Since $D$ has degree $9$, its coefficient vanishes for $n\ge10$; with $f_m=0$ for $m\lt0$ this leaves a homogeneous recurrence of order $18$, valid from $n=10$ on. The numerator matters too: a denominator of degree $d$ alone would not bound the order from the start of a sequence. Written out,

$$f_n=11\sum_{d=1}^{9}f_{n-d}-10\sum_{i,j=1}^{9}f_{n-i-j},$$

and we verified this for $10\le n\le60$. The OEIS entry for $f$ records the same order $18$. ^[4]

6. A closed form when digits may exceed $9$

Drop the bound on the parts, giving a composition $(c_1,\ldots,c_k)$ the value $\sum_j c_j10^{k-j}$ as if each $c_j$ were a digit. Then $P$ and $D$ become $x/(1-x)$ and $x/(1-x)^2$, the factors $(1-x)^2$ cancel, and the quotient collapses to $x/\bigl((1-2x)(1-11x)\bigr)$, whose coefficients are $(11^n-2^n)/9$. ^[5] No composition of an $n\le9$ has a part above $9$, so

$$f_n=\frac{11^n-2^n}{9}\qquad(1\le n\le9),$$

for instance $f_5=(161051-32)/9=17891$. From $n=10$ the two part company: the formula gives $2881935953$ against $f_{10}=2881935943$, a difference of $10$ — exactly the value of the banned composition $(10)$. The bounded problem has no such two-term formula, so we evaluate its order-$18$ recurrence by matrix powers instead.

7. Why the recurrence cannot simply be iterated

The largest index we need is

$$13^{17}=8\,650\,415\,919\,381\,337\,933,$$

a $63$-bit number, so stepping once per index would take more than $10^{19}$ additions. Nor can the integers themselves be touched: a zeroless integer of digit sum $13^{17}$ has at least $\lceil 13^{17}/9\rceil=961\,157\,324\,375\,704\,215$ digits, about $9.6\cdot10^{17}$ of them. Only the algebra of the recurrence is usable.

8. One matrix for both recurrences

Both recurrences look back exactly nine steps, so the nine most recent values of each sequence are a complete state. Collect them in one column vector of length $18$,

$$v_n=(a_n,\ldots,a_{n-8},f_n,\ldots,f_{n-8})^{\mathrm T},$$

and read the two recurrences as instructions for the next state:

$$a_{n+1}=\sum_{i=0}^{8}a_{n-i},$$

$$f_{n+1}=\sum_{i=0}^{8}\bigl((i+1)a_{n-i}+10f_{n-i}\bigr).$$

The remaining sixteen entries of $v_{n+1}$ are entries of $v_n$ shifted by one place, so there is a single matrix $T$ with

$$v_{n+1}=Tv_n,\qquad v_n=T^nv_0.$$

Its first row is nine ones; rows $2$ to $9$ each hold a single $1$ shifting the count window; row $10$ holds $1,2,\ldots,9$ on the count half and nine copies of $10$ on the sum half; rows $11$ to $18$ shift the sum window. Only $43$ of the $324$ entries are nonzero. The starting state is $v_0=(1,0,\ldots,0)^{\mathrm T}$, from $a_0=1$ and $f_0=0$, and $f_n$ is the tenth entry of $v_n$. This is the companion-matrix picture of a constant-recursive sequence, blocked by the one-way coupling of Section 3. ^[6]

9. Raising the matrix by repeated squaring

Computing $T^n$ by $n-1$ multiplications would be no better than iterating. Instead square repeatedly, $T^{2^0},T^{2^1},T^{2^2},\ldots$, and combine those powers whose bit is set in $n$: $\lfloor\log_2 n\rfloor$ squarings plus one multiplication per further one-bit. The method works for matrices and modular arithmetic alike, ^[7] and powers of one matrix commute, so the order of the factors is irrelevant. Take the first target, $13=1+4+8$: build $T,T^2,T^4,T^8$, then apply $T$, $T^4$ and $T^8$ to $v_0$ in turn to obtain $T^{13}v_0$ — three vector products, and $T^{13}$ itself is never formed.

One refinement matters: multiplying two $18\times18$ matrices costs $18^3$ multiply-adds, while applying one to a vector costs $18^2$. As every query starts from the same $v_0$, the squarings are done once and each query only applies matrices to a vector. The seventeen exponents $13^i$ carry between $3$ and $32$ one-bits, $275$ altogether, so $275$ vector products answer all of them.

10. Working modulo $10^9$, and what it does to the matrix

Only the last nine digits of $S$ are wanted, and every operation in sight is an addition or a multiplication by a constant. Reduction modulo $10^9$ commutes with both, so every intermediate value may be reduced immediately. No division or modular inverse occurs anywhere, which is why it does not matter that $10^9=2^9\cdot5^9$ is far from prime.

The modulus also makes the block that carries old sums to new sums vanish, which is why the sparsity of $T$ survives squaring. Let $C_t(g)$ count the compositions of $g$ into exactly $t$ parts from $\{1,\ldots,9\}$, with $C_0(0)=1$, and set

$$c_g=\sum_{t\ge0}C_t(g)\,10^{t}.$$

Appending a suffix of $t$ digits multiplies a prefix by $10^{t}$, so $c_g$ is the coefficient of the newest sum $f_n$ in $f_{n+g}$; the values $c_1=10$, $c_2=110$, $c_3=1210$ match the matrix entries, checked for $g\le29$. Older entries need a restriction, since the nine sums in the state are independent starting values: a chain of steps from $f_{n-i}$ up to $f_{n+k-j}$ must not pass through another of them, so its first step must exceed $i$. The coefficient from $f_{n-1}$ to $f_{n+1}$ is thus $10$, not $c_2$.

Each allowed chain of total rise $g$ still takes at least $\lceil g/9\rceil$ steps, every step carrying a factor $10$, so each entry is divisible by $10^{\lceil g/9\rceil}$ and vanishes modulo $10^9$ once $g\ge73$. Here $g=k-j+i\ge k-8$, so the whole block vanishes for $k\ge81$. At $k=80$ only $(j,i)=(8,0)$ leaves $g=72$, whose one eight-step chain is eight nines and contributes $10^{8}$: exactly one surviving entry, as we confirmed. The binary powers therefore fall from $243$ nonzero entries at $T^{16}$ and $T^{32}$ to $162$ for every $T^{2^b}$ with $b\ge7$, each row then holding at most nine.

How the Code Works

The C++ program ^[8] transcribes Sections 8 to 10 directly. build_transition writes $T$ row by row as described there; multiply is the schoolbook $18^3$ product, reducing modulo kMod $=10^9$ after every term, skipping a term as soon as either factor is zero — the skips are what the sparsity of Section 10 pays for — and forming each product in __int128 first. ^[9] precompute_powers squares repeatedly to store $T^{2^b}$ for $b=0,\ldots,63$.

f_of_n sets the state to $v_0$, walks the bits of $n$ upward from the least significant, calls apply_matrix_vector once per set bit and returns state[9], which is $f_n$ modulo $10^9$; main multiplies a running p13 by $13$ seventeen times and adds the residues. Two details are slack rather than wrong: identity_matrix is never called, and $T^{2^{63}}$ is built although $13^{17}$ occupies $63$ bits, so index $62$ is the highest ever used.

run_checkpoints is the self-test: unless --skip-checkpoints is given, it verifies $f(5)=17891$, then builds the count and sum arrays by the recurrences of Sections 2 and 3 for $1\le n\le29$ and compares each f[n] with f_of_n — the counts feed that computation but are not themselves compared.

The three files differ in how much they factor out. Python keeps build_transition, multiply, apply_matrix_vector and f_of_n, squares inside solve, and has no identity-matrix helper and no checkpoints; Java keeps only buildTransition, multiply and applyMatrixVector, holding the squaring and every bit loop in solve. Python's integers are unbounded, ^[9] so nothing there can overflow. Java uses long[][]: in multiply an accumulator stays below $(10^9-1)^2+(10^9-1)$, well inside the signed $64$-bit range, but applyMatrixVector adds up to eighteen products of two residues before reducing, and $18(10^9-1)^2$ exceeds $2^{63}-1$; the next section reports what we measured.

Complexity and Verification

With state dimension $d=18$ the cost is $63$ squarings at $d^3=5832$ multiply-adds each plus $d^2=324$ per set bit, that is $O(d^3\log n)$ for the largest exponent $n$. Densely that would be $63\cdot5832+275\cdot324=456\,516$ multiply-adds; with the zero skips we measured $93\,714$ for the squarings and $40\,713$ for the queries, about $134$ thousand in all, and the C++ checkpoints add $73$ more matrix-vector applications plus a small direct recurrence. The stored C++ entries occupy $165\,888$ bytes, that is $162$ KiB, excluding container overhead; Java holds the same long payload, Python references to arbitrary-precision integers and so more. All three finish in a small fraction of a second. Since $f$ obeys one order-$18$ recurrence (Section 5), reducing $x^n$ modulo its characteristic polynomial would cost $O(d^2\log n)$; at $d=18$ that is immaterial, but we implemented it as an independent check.

The following were verified for this article.

The recurrences of Sections 2 and 3 agree with a direct enumeration of every zeroless integer of digit sum $n$ for $1\le n\le15$, all $16\,272$ of them at $n=15$.
The nine values $f_1,\ldots,f_9$ match the closed form $(11^n-2^n)/9$ of Section 6, the statement's checkpoint $f_5=17891$ among them.
$f_{20}=74\,749\,998\,860\,563\,784\,655$ agrees with the published term of the OEIS entry for $f$. ^[4]
The matrix method equals the plain dynamic programme modulo $10^9$ for every $1\le n\le200$, and the order-$18$ recurrence of Section 5 holds for $10\le n\le60$.
An independent program that powers $x^n$ modulo the characteristic polynomial agrees with the matrix for $0\le n\lt400$ and at all seventeen exponents $13^i$, and reaches the same total. It is not part of the published solutions.
The C++, Python and Java programs print the same value, and the C++ self-test passes.

The Java accumulator deserved a closer look, so we summed the rows of all $64$ binary powers. For three of them — $T^{16}$, $T^{32}$ and $T^{64}$, all used here — the largest row sum times $10^9-1$ exceeds $2^{63}-1$, so a state whose entries were all $10^9-1$ would overflow applyMatrixVector; the largest such hypothetical accumulator is $1.13\cdot10^{19}$. On the states that actually arise the largest accumulator is $5.38\cdot10^{18}$, a factor $1.71$ below the limit, so the Java port is correct — but on these particular numbers, not by a bound, and Java's integer operators signal no overflow. ^[9] From $T^{128}$ on the sum-to-sum block has vanished (Section 10), every row holds at most nine entries, and $9(10^9-1)^2$ stays below $2^{63}-1$, so those powers are provably safe. C++ accumulates in $128$ bits, Python without bound.

Footnotes and References

The references below supply the background facts used above. The generating function for f, the order-18 recurrence, the valuation of the sum-to-sum coefficients and the operation counts are derived in this article.

Project Euler 377 — Sum of Digits - Experience #13. The official statement lists the 16 positive integers with no zero digit and digital sum 5, records that they add up to 17891, defines f(n) as the sum of all positive integers without a zero digit whose digital sum is n, and asks for the last 9 digits of the sum of f(13^i) over i = 1 to 17.
Wikipedia — Composition (combinatorics). A composition writes n as an ordered sum of positive integers; each n has 2^(n−1) of them, one per choice of cutting or joining each of the n−1 gaps; and for parts drawn from a set A, the k-part compositions are counted by the coefficient of x^n in the k-th power of the sum of x^a over a in A. Sections 1 and 5 use all three.
OEIS A258800 — zeroless integers of digital sum n. "The number of zeroless decimal numbers whose digital sum is n": terms 1, 2, 4, 8, 16, 32, 64, 128, 256, 511, 1021, 2040, 4076 from n = 1, and a generating function of denominator 1 − x − x² − … − x⁹. Its formula a(n) = A104144(n + 8) identifies it with the nine-step Fibonacci numbers. These are the counts of Section 2, save that the entry sets a(0) = 0 where our empty-string convention sets 1.
OEIS A211072 — sums of zeroless integers of digital sum n. "Sum of numbers with no 0 decimal digits whose sum of digits equals n": 1, 13, 147, 1625, 17891, … and 74749998860563784655 at n = 20, with a rational generating function, indexed under linear recurrences of order 18, and noted as different from A016135. It is the sequence f, confirming Sections 5 and 6.
OEIS A016135 — expansion of 1/((1 − 2x)(1 − 11x)). Its closed form is a(n) = (11^(n+1) − 2^(n+1))/9, with terms 1, 13, 147, 1625, 17891, 196833, 2165227, 23817625, 261994131, 2881935953. The first nine are f(1) to f(9); the tenth is where the unrestricted model of Section 6 overshoots f(10) = 2881935943 by 10.
Wikipedia — Constant-recursive sequence. A sequence satisfies a homogeneous linear recurrence with constant coefficients of order at most d, valid from n = d on, exactly when its generating function can be written p(x)/q(x) with q(0) = 1, the degree of p below d and the degree of q at most d; such a sequence also has the form u·Aⁿ·v for a d×d companion matrix A. Sections 5 and 8 use both connections.
Wikipedia — Exponentiation by squaring. Reading the exponent in binary needs ⌊log₂ n⌋ squarings and one further multiplication per one-bit beyond the first; the method applies in any semigroup, hence to matrices, and is presented there as the standard way to take large powers modulo a number. This is Section 9.
C++ — Euler377.cpp. The linked immutable revision contains build_transition, multiply, apply_matrix_vector, precompute_powers, f_of_n and run_checkpoints, with kSize 18 and kMod 1000000000. It is the implementation record for this article; the polynomial cross-check above is extra validation, not a test in that file.
GCC manual — 128-bit integers documents __int128 as an extension on targets with a 128-bit integer mode; it holds the C++ products before reduction. The Java Language Specification, §4.2 gives long the range −9223372036854775808 to 9223372036854775807 and states that the integer operators do not indicate overflow or underflow in any way, so an accumulator past the limit wraps silently. Python — numeric types states that integers have unlimited precision.

Problem 377 source code

C++

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

namespace {

using i64 = std::int64_t;
using u64 = std::uint64_t;

constexpr i64 kMod = 1000000000LL;
constexpr int kSize = 18;

using Matrix = std::array<std::array<i64, kSize>, kSize>;
using Vector = std::array<i64, kSize>;

Matrix zero_matrix() {
    Matrix m{};
    for (int i = 0; i < kSize; ++i) {
        for (int j = 0; j < kSize; ++j) {
            m[i][j] = 0;
        }
    }
    return m;
}

Matrix identity_matrix() {
    Matrix m = zero_matrix();
    for (int i = 0; i < kSize; ++i) {
        m[i][i] = 1;
    }
    return m;
}

Matrix multiply(const Matrix& a, const Matrix& b) {
    Matrix c = zero_matrix();
    for (int i = 0; i < kSize; ++i) {
        for (int k = 0; k < kSize; ++k) {
            const i64 aik = a[i][k];
            if (aik == 0) {
                continue;
            }
            for (int j = 0; j < kSize; ++j) {
                if (b[k][j] == 0) {
                    continue;
                }
                c[i][j] = (c[i][j] + static_cast<__int128>(aik) * b[k][j]) % kMod;
            }
        }
    }
    return c;
}

Vector apply_matrix_vector(const Matrix& a, const Vector& v) {
    Vector out{};
    for (int i = 0; i < kSize; ++i) {
        __int128 sum = 0;
        for (int j = 0; j < kSize; ++j) {
            if (a[i][j] == 0 || v[j] == 0) {
                continue;
            }
            sum += static_cast<__int128>(a[i][j]) * v[j];
        }
        out[i] = static_cast<i64>(sum % kMod);
    }
    return out;
}

Matrix build_transition() {
    Matrix t = zero_matrix();

    // a_{n+1} = a_n + a_{n-1} + ... + a_{n-8}
    for (int i = 0; i < 9; ++i) {
        t[0][i] = 1;
    }
    for (int i = 1; i < 9; ++i) {
        t[i][i - 1] = 1;
    }

    // f_{n+1} = 10*(f_n+...+f_{n-8}) + 1*a_n + 2*a_{n-1} + ... + 9*a_{n-8}
    for (int i = 0; i < 9; ++i) {
        t[9][i] = i + 1;
        t[9][9 + i] = 10;
    }
    for (int i = 10; i < 18; ++i) {
        t[i][i - 1] = 1;
    }

    return t;
}

std::vector<Matrix> precompute_powers(const Matrix& base) {
    std::vector<Matrix> powers;
    powers.reserve(64);
    powers.push_back(base);
    for (int i = 1; i < 64; ++i) {
        powers.push_back(multiply(powers.back(), powers.back()));
    }
    return powers;
}

i64 f_of_n(const u64 n, const std::vector<Matrix>& powers) {
    // state at n=0:
    // [a_0, a_{-1},...,a_{-8}, f_0, f_{-1},...,f_{-8}] with a_0=1, others 0.
    Vector state{};
    for (int i = 0; i < kSize; ++i) {
        state[i] = 0;
    }
    state[0] = 1;

    u64 e = n;
    int bit = 0;
    while (e > 0) {
        if ((e & 1ULL) != 0ULL) {
            state = apply_matrix_vector(powers[static_cast<std::size_t>(bit)], state);
        }
        e >>= 1ULL;
        ++bit;
    }

    // f_n is at index 9.
    return state[9];
}

bool run_checkpoints(const std::vector<Matrix>& powers) {
    if (f_of_n(5ULL, powers) != 17891LL) {
        std::cerr << "Checkpoint failed: f(5)\n";
        return false;
    }

    // Cross-check with direct DP for small n.
    std::vector<i64> a(30, 0);
    std::vector<i64> f(30, 0);
    a[0] = 1;
    for (int n = 1; n < 30; ++n) {
        for (int d = 1; d <= 9; ++d) {
            if (n - d < 0) {
                continue;
            }
            a[n] = (a[n] + a[n - d]) % kMod;
            f[n] = (f[n] + 10LL * f[n - d] + static_cast<i64>(d) * a[n - d]) % kMod;
        }
    }
    for (int n = 1; n < 30; ++n) {
        if (f_of_n(static_cast<u64>(n), powers) != f[n]) {
            std::cerr << "Checkpoint failed: matrix/DP mismatch at n=" << n << '\n';
            return false;
        }
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    bool skip_checkpoints = false;
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            skip_checkpoints = true;
        } else {
            std::cerr << "Unknown argument: " << arg << '\n';
            return 1;
        }
    }

    const Matrix transition = build_transition();
    const std::vector<Matrix> powers = precompute_powers(transition);

    if (!skip_checkpoints && !run_checkpoints(powers)) {
        return 2;
    }

    i64 answer = 0;
    u64 p13 = 1ULL;
    for (int i = 1; i <= 17; ++i) {
        p13 *= 13ULL;
        answer += f_of_n(p13, powers);
        answer %= kMod;
    }

    std::cout << answer << '\n';
    return 0;
}

Python

def solve():
    kMod = 1000000000
    kSize = 18

    def zero_matrix():
        return [[0]*kSize for _ in range(kSize)]

    def multiply(a, b):
        c = zero_matrix()
        for i in range(kSize):
            for k in range(kSize):
                if a[i][k] == 0: continue
                for j in range(kSize):
                    if b[k][j] == 0: continue
                    c[i][j] = (c[i][j] + a[i][k] * b[k][j]) % kMod
        return c

    def apply_matrix_vector(a, v):
        out = [0]*kSize
        for i in range(kSize):
            s = 0
            for j in range(kSize):
                if a[i][j] == 0 or v[j] == 0: continue
                s += a[i][j] * v[j]
            out[i] = s % kMod
        return out

    def build_transition():
        t = zero_matrix()
        for i in range(9):
            t[0][i] = 1
        for i in range(1, 9):
            t[i][i - 1] = 1
        for i in range(9):
            t[9][i] = i + 1
            t[9][9 + i] = 10
        for i in range(10, 18):
            t[i][i - 1] = 1
        return t

    transition = build_transition()
    powers = [transition]
    for _ in range(1, 64):
        powers.append(multiply(powers[-1], powers[-1]))

    def f_of_n(n):
        state = [0]*kSize
        state[0] = 1
        e = n
        bit = 0
        while e > 0:
            if (e & 1) != 0:
                state = apply_matrix_vector(powers[bit], state)
            e >>= 1
            bit += 1
        return state[9]

    ans = 0
    p13 = 1
    for _ in range(1, 18):
        p13 *= 13
        ans = (ans + f_of_n(p13)) % kMod

    return str(ans)

if __name__ == '__main__':
    print(solve())

Java

public class Euler377 {
    private static final long kMod = 1000000000L;
    private static final int kSize = 18;

    private static long[][] zeroMatrix() {
        return new long[kSize][kSize];
    }

    private static long[][] multiply(long[][] a, long[][] b) {
        long[][] c = zeroMatrix();
        for (int i = 0; i < kSize; ++i) {
            for (int k = 0; k < kSize; ++k) {
                if (a[i][k] == 0)
                    continue;
                for (int j = 0; j < kSize; ++j) {
                    if (b[k][j] == 0)
                        continue;
                    c[i][j] = (c[i][j] + a[i][k] * b[k][j]) % kMod;
                }
            }
        }
        return c;
    }

    private static long[] applyMatrixVector(long[][] a, long[] v) {
        long[] out = new long[kSize];
        for (int i = 0; i < kSize; ++i) {
            long sum = 0;
            for (int j = 0; j < kSize; ++j) {
                if (a[i][j] == 0 || v[j] == 0)
                    continue;
                sum += a[i][j] * v[j];
            }
            out[i] = sum % kMod;
        }
        return out;
    }

    private static long[][] buildTransition() {
        long[][] t = zeroMatrix();
        for (int i = 0; i < 9; ++i) {
            t[0][i] = 1;
        }
        for (int i = 1; i < 9; ++i) {
            t[i][i - 1] = 1;
        }
        for (int i = 0; i < 9; ++i) {
            t[9][i] = i + 1;
            t[9][9 + i] = 10;
        }
        for (int i = 10; i < 18; ++i) {
            t[i][i - 1] = 1;
        }
        return t;
    }

    public static String solve() {
        long[][] transition = buildTransition();
        long[][][] powers = new long[64][][];
        powers[0] = transition;
        for (int i = 1; i < 64; ++i) {
            powers[i] = multiply(powers[i - 1], powers[i - 1]);
        }

        long ans = 0;
        long p13 = 1;
        for (int i = 1; i <= 17; ++i) {
            p13 *= 13L;

            long[] state = new long[kSize];
            state[0] = 1;
            long e = p13;
            int bit = 0;
            while (e > 0) {
                if ((e & 1L) != 0L) {
                    state = applyMatrixVector(powers[bit], state);
                }
                e >>= 1L;
                ++bit;
            }

            ans = (ans + state[9]) % kMod;
        }
        return String.valueOf(ans);
    }

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