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

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)$.

Adding the contributions of the nine possible last digits gives the coupled recurrence

$$F(n)=\sum_{d=1}^{9}\left(10F(n-d)+d\,A(n-d)\right),$$

with initial conditions

$$F(0)=0,\qquad F(n)=0 \text{ for } n\lt 0.$$

This pair of recurrences is exactly what the C++, Python, and Java solutions implement.

Step 3: Small example at $n=5$

The valid numbers are $5$; $14,23,32,41$; $113,122,131,212,221,311$; $1112,1121,1211,2111$; and $11111$. There are $16$ such numbers, so $A(5)=16$, and their total is

$$F(5)=17891.$$

The C++ program uses this identity as a checkpoint before it starts the large matrix-power computations.

Step 4: Turn the recurrence into an 18-dimensional linear system

Both recurrences only look back nine steps, so one state vector stores everything needed for the next update:

$$v_n=\bigl(A(n),A(n-1),\dots,A(n-8),F(n),F(n-1),\dots,F(n-8)\bigr)^T.$$

From the formulas above we obtain

$$A(n+1)=A(n)+A(n-1)+\cdots+A(n-8),$$

$$F(n+1)=\sum_{i=0}^{8}(i+1)A(n-i)+10\sum_{i=0}^{8}F(n-i).$$

All other coordinates of $v_{n+1}$ are simple shifts of the previous window. Hence there exists a fixed $18\times18$ matrix $T$ such that

$$v_{n+1}=T\,v_n,\qquad v_n=T^n v_0,$$

with initial state

$$v_0=(1,0,\dots,0,0,\dots,0)^T.$$

The required value $F(n)$ is the 10th component of $v_n$, which is why the implementations read state index $9$.

Step 5: Use binary exponentiation for the indices $13^k$

The target indices are $13,13^2,\dots,13^{17}$, so iterating the recurrence step by step is still infeasible. Instead the code precomputes

$$T^{2^0},T^{2^1},T^{2^2},\dots$$

by repeated squaring modulo $10^9$. For each exponent $n$, the binary expansion of $n$ tells us which precomputed powers must be applied to $v_0$. Because the recurrence is linear, reducing every intermediate result modulo $10^9$ is mathematically valid.

How the Code Works

build_transition() constructs the sparse $18\times18$ matrix described above. The first row contains nine ones for the count recurrence, the next eight rows shift the count window, the 10th row contains coefficients $1,2,\dots,9$ for the count block and nine copies of $10$ for the value-sum block, and the last eight rows shift the value-sum window.

The C++ version precomputes 64 matrix powers because $13^{17}$ fits within 64 bits, then f_of_n() applies only the powers corresponding to set bits of $n$. The Python and Java versions use the same transition and the same binary-exponentiation idea.

Complexity Analysis

Let $m$ be the largest queried index. Precomputing the matrix powers costs $O(18^3\log m)$ time and $O(18^2\log m)$ memory. After that, one evaluation of $F(n)$ costs only $O(18^2\log n)$ time because the query phase multiplies matrices by a vector, not by another matrix. For this problem there are only 17 queries, so the total runtime is dominated by a small constant number of $18\times18$ operations.

Footnotes and References

Problem page: https://projecteuler.net/problem=377
Matrix exponentiation for linear recurrences: cp-algorithms
Generating functions and coefficient extraction: Wikipedia — Generating function
Digit sums and compositions into parts $1$ through $9$: Wikipedia — Composition

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