Problem 778: Freshman's Product

View on Project Euler

Project Euler Problem 778 Solution

EulerSolve provides an optimized solution for Project Euler Problem 778, Freshman's Product, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Write every nonnegative integer in decimal, padding with leading zeros when necessary so that all positions line up. The freshman's product multiplies matching digits and keeps only the last digit at each position, with no carries between positions. If $$d_j(n)=\left\lfloor\frac{n}{10^j}\right\rfloor \bmod 10,$$ then for two numbers $$a\odot b=\sum_{j\ge 0}\left(d_j(a)d_j(b)\bmod 10\right)10^j.$$ Because multiplication modulo \(10\) is associative at each digit, the same definition extends to \(R\) factors: $$n_1\odot n_2\odot \cdots \odot n_R=\sum_{j\ge 0}\left(\prod_{t=1}^R d_j(n_t)\bmod 10\right)10^j.$$ The quantity to compute is $$F(R,M)=\sum_{0\le n_1,\dots,n_R\le M} n_1\odot n_2\odot \cdots \odot n_R \pmod{10^9+9}.$$ A direct enumeration would require \((M+1)^R\) ordered tuples, so the solution instead counts the contribution of each decimal position separately and compresses the repeated digit multiplication into a fixed \(10\times 10\) transition matrix. Mathematical Approach The crucial observation is that the freshman's product never creates carries. Each decimal place can therefore be analyzed independently and then recombined with its place value. Step 1: Separate the Decimal Positions Fix one place \(p=10^j\)....

Detailed mathematical approach

Problem Summary

Write every nonnegative integer in decimal, padding with leading zeros when necessary so that all positions line up. The freshman's product multiplies matching digits and keeps only the last digit at each position, with no carries between positions. If

$$d_j(n)=\left\lfloor\frac{n}{10^j}\right\rfloor \bmod 10,$$

then for two numbers

$$a\odot b=\sum_{j\ge 0}\left(d_j(a)d_j(b)\bmod 10\right)10^j.$$

Because multiplication modulo \(10\) is associative at each digit, the same definition extends to \(R\) factors:

$$n_1\odot n_2\odot \cdots \odot n_R=\sum_{j\ge 0}\left(\prod_{t=1}^R d_j(n_t)\bmod 10\right)10^j.$$

The quantity to compute is

$$F(R,M)=\sum_{0\le n_1,\dots,n_R\le M} n_1\odot n_2\odot \cdots \odot n_R \pmod{10^9+9}.$$

A direct enumeration would require \((M+1)^R\) ordered tuples, so the solution instead counts the contribution of each decimal position separately and compresses the repeated digit multiplication into a fixed \(10\times 10\) transition matrix.

Mathematical Approach

The crucial observation is that the freshman's product never creates carries. Each decimal place can therefore be analyzed independently and then recombined with its place value.

Step 1: Separate the Decimal Positions

Fix one place \(p=10^j\). For an ordered \(R\)-tuple \((n_1,\dots,n_R)\), the digit contributed at that place is

$$g_p(n_1,\dots,n_R)=\left(\prod_{t=1}^R d_j(n_t)\right)\bmod 10.$$

So the whole sum splits as

$$F(R,M)=\sum_{\substack{j\ge 0\\10^j\le M}} 10^j S_j \pmod{10^9+9},$$

where

$$S_j=\sum_{0\le n_1,\dots,n_R\le M} g_{10^j}(n_1,\dots,n_R).$$

All higher positions vanish automatically, because every number in \([0,M]\) has digit \(0\) there.

Step 2: Count How Often Each Digit Appears at One Place

For a fixed place \(p=10^j\), let \(c_d(p)\) be the number of integers in \([0,M]\) whose digit at that place equals \(d\). Set

$$N=M+1,\qquad q=\left\lfloor\frac{N}{10p}\right\rfloor,\qquad u=N\bmod(10p).$$

Each complete block of length \(10p\) contributes exactly \(p\) copies of every digit, so the \(q\) full blocks contribute \(qp\). The remaining tail contributes a partial block, giving

$$c_d(p)=qp+\min\bigl(\max(u-dp,0),\,p\bigr),\qquad 0\le d\le 9.$$

This formula automatically handles leading zeros. For example, if \(M=27\) and \(p=10\), then the tens digit counts are

$$c_0(10)=10,\qquad c_1(10)=10,\qquad c_2(10)=8,\qquad c_d(10)=0\ \text{for }d\ge 3,$$

because the numbers \(0\) to \(9\) contribute a leading zero in the tens place.

Step 3: Model Repeated Multiplication with 10 Residue States

At one fixed decimal place, only the current product modulo \(10\) matters. Let \(v_t^{(j)}(r)\) be the number of ordered choices of \(t\) numbers from \([0,M]\) whose digitwise product at place \(10^j\) equals \(r\in\{0,\dots,9\}\).

Before choosing any number, the product is the multiplicative identity, so

$$v_0^{(j)}(1)=1,\qquad v_0^{(j)}(r)=0\ \text{for }r\ne 1.$$

If the current residue is \(r\) and the next chosen digit is \(d\), the new residue is

$$r'\equiv rd \pmod{10}.$$

This gives a linear transition matrix \(T_j\) of size \(10\times 10\):

$$T_j(r',r)=\sum_{\substack{0\le d\le 9\\ rd\equiv r' \pmod{10}}} c_d(10^j).$$

Then the state update is simply

$$v_{t+1}^{(j)}=T_j\,v_t^{(j)}.$$

These entries are counts, not probabilities: every column sums to \(M+1\), because from any current residue there are \(M+1\) possible next numbers.

Step 4: Jump Directly to \(R\) Factors by Matrix Exponentiation

Let \(e_1\) be the basis vector concentrated at residue \(1\). After choosing \(R\) numbers,

$$v_R^{(j)}=T_j^R e_1.$$

The total digit contribution at place \(10^j\) is therefore

$$S_j=\sum_{r=0}^9 r\,v_R^{(j)}(r)=\sum_{r=0}^9 r\,(T_j^R e_1)_r.$$

Since \(R\) can be very large, binary exponentiation computes \(T_j^R\) in \(O(10^3\log R)\) operations for that place.

Step 5: Recombine All Decimal Places

Putting everything together gives the exact formula implemented by the program:

$$\boxed{F(R,M)=\sum_{\substack{j\ge 0\\10^j\le M}} 10^j \sum_{r=0}^9 r\,(T_j^R e_1)_r \pmod{10^9+9}.}$$

The only place-dependent ingredient is the digit-count vector \(c_d(10^j)\); once those ten counts are known, the rest is the same \(10\)-state linear algebra for every decimal position.

Worked Example: \((R,M)=(2,7)\)

Here only the units place exists, because every number from \(0\) to \(7\) has just one relevant digit. The digit counts are

$$c_0(1)=c_1(1)=\cdots=c_7(1)=1,\qquad c_8(1)=c_9(1)=0.$$

So the desired sum is simply

$$F(2,7)=\sum_{a=0}^7\sum_{b=0}^7 (ab\bmod 10).$$

Evaluating row by row gives

$$\begin{aligned} a=0&:&&0,\\ a=1&:&&0+1+2+3+4+5+6+7=28,\\ a=2&:&&0+2+4+6+8+0+2+4=26,\\ a=3&:&&0+3+6+9+2+5+8+1=34,\\ a=4&:&&0+4+8+2+6+0+4+8=32,\\ a=5&:&&0+5+0+5+0+5+0+5=20,\\ a=6&:&&0+6+2+8+4+0+6+2=28,\\ a=7&:&&0+7+4+1+8+5+2+9=36. \end{aligned}$$

Hence

$$0+28+26+34+32+20+28+36=204,$$

so \(F(2,7)=204\), matching the implementation checkpoint.

How the Code Works

The C++, Python, and Java implementations follow the same pipeline. They iterate over the place values \(1,10,100,\dots\) up to the largest nonzero digit of \(M\). For each place they compute the ten digit frequencies in \([0,M]\) from the block formula above, assemble the \(10\times 10\) transition matrix on residues modulo \(10\), raise that matrix to the \(R\)-th power by binary exponentiation, and read the weighted digit sum from the column corresponding to the initial residue \(1\). That place contribution is then multiplied by the place value and added to the running answer modulo \(10^9+9\). The three implementations differ only in language syntax and integer handling; the mathematics and the sequence of steps are identical.

Complexity Analysis

Let \(L=\lfloor\log_{10} M\rfloor+1\) be the number of decimal positions that must be processed. For one position, computing the digit counts costs \(O(10)\), building the transition matrix costs \(O(10^2)\), and exponentiating a \(10\times 10\) matrix costs \(O(10^3\log R)\). Therefore the total running time is

$$O\bigl(L\cdot 10^3\log R\bigr),$$

which is effectively \(O(L\log R)\) because the state space size \(10\) is fixed. The memory usage is \(O(10^2)\) for a few small matrices and digit-count vectors.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=778
  2. Positional notation: Wikipedia — Positional notation
  3. Modular arithmetic: Wikipedia — Modular arithmetic
  4. Exponentiation by squaring: Wikipedia — Exponentiation by squaring
  5. Matrix multiplication: Wikipedia — Matrix multiplication

Problem 778 source code

C++

#include <array>
#include <cassert>
#include <cstdint>
#include <iostream>

namespace {

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

constexpr i64 MOD = 1'000'000'009LL;
constexpr int D = 10;

using Matrix = std::array<std::array<i64, D>, D>;

Matrix mat_mul(const Matrix& A, const Matrix& B) {
    Matrix C{};
    for (int i = 0; i < D; ++i) {
        for (int k = 0; k < D; ++k) {
            const i64 aik = A[i][k];
            if (aik == 0) {
                continue;
            }
            for (int j = 0; j < D; ++j) {
                C[i][j] = static_cast<i64>((C[i][j] + static_cast<__int128>(aik) * B[k][j]) % MOD);
            }
        }
    }
    return C;
}

Matrix mat_pow(Matrix base, u64 exp) {
    Matrix result{};
    for (int i = 0; i < D; ++i) {
        result[i][i] = 1;
    }

    while (exp > 0) {
        if (exp & 1ULL) {
            result = mat_mul(base, result);
        }
        base = mat_mul(base, base);
        exp >>= 1ULL;
    }
    return result;
}

std::array<i64, D> digit_counts(u64 M, u64 p10) {
    const u64 N = M + 1;
    const u64 block = 10ULL * p10;
    const u64 full = N / block;
    const u64 rem = N % block;

    std::array<i64, D> c{};
    for (int d = 0; d < D; ++d) {
        u64 cnt = full * p10;
        const u64 shift = static_cast<u64>(d) * p10;
        if (rem > shift) {
            const u64 extra = rem - shift;
            cnt += (extra < p10) ? extra : p10;
        }
        c[d] = static_cast<i64>(cnt % MOD);
    }
    return c;
}

i64 contribution_for_position(u64 R, u64 M, u64 p10) {
    const auto c = digit_counts(M, p10);

    Matrix T{};
    for (int from = 0; from < D; ++from) {
        for (int d = 0; d < D; ++d) {
            const int to = (from * d) % 10;
            T[to][from] += c[d];
            if (T[to][from] >= MOD) {
                T[to][from] -= MOD;
            }
        }
    }

    const Matrix P = mat_pow(T, R);

    i64 s = 0;
    for (int r = 0; r < D; ++r) {
        s = static_cast<i64>((s + static_cast<__int128>(r) * P[r][1]) % MOD);
    }
    return s;
}

i64 F_mod(u64 R, u64 M) {
    i64 ans = 0;
    i64 place_mod = 1;

    for (u64 p10 = 1; p10 <= M; p10 *= 10ULL) {
        const i64 s = contribution_for_position(R, M, p10);
        ans = static_cast<i64>((ans + static_cast<__int128>(place_mod) * s) % MOD);
        place_mod = (place_mod * 10LL) % MOD;
    }

    return ans;
}

}  // namespace

int main() {
    assert(F_mod(2, 7) == 204);
    assert(F_mod(23, 76) == 5'870'548);

    std::cout << F_mod(234'567, 765'432) << '\n';
    return 0;
}

Python

MOD = 1000000009
D = 10

def mat_mul(A, B):
    C = [[0] * D for _ in range(D)]
    for i in range(D):
        for k in range(D):
            aik = A[i][k]
            if aik == 0: continue
            for j in range(D):
                C[i][j] = (C[i][j] + aik * B[k][j]) % MOD
    return C

def mat_pow(base, exp):
    result = [[1 if i == j else 0 for j in range(D)] for i in range(D)]
    while exp > 0:
        if exp & 1:
            result = mat_mul(base, result)
        base = mat_mul(base, base)
        exp >>= 1
    return result

def digit_counts(M, p10):
    N = M + 1
    block = 10 * p10
    full = N // block
    rem = N % block
    
    c = [0] * D
    for d in range(D):
        cnt = full * p10
        shift = d * p10
        if rem > shift:
            extra = rem - shift
            cnt += extra if extra < p10 else p10
        c[d] = cnt % MOD
    return c

def contribution_for_position(R, M, p10):
    c = digit_counts(M, p10)
    T = [[0] * D for _ in range(D)]
    for frm in range(D):
        for d in range(D):
            to = (frm * d) % 10
            T[to][frm] = (T[to][frm] + c[d]) % MOD
            
    P = mat_pow(T, R)
    s = 0
    for r in range(D):
        s = (s + r * P[r][1]) % MOD
    return s

def f_mod(R, M):
    ans = 0
    place_mod = 1
    p10 = 1
    while p10 <= M:
        s = contribution_for_position(R, M, p10)
        ans = (ans + place_mod * s) % MOD
        place_mod = (place_mod * 10) % MOD
        p10 *= 10
    return ans

def solve():
    return str(f_mod(234567, 765432))

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

Java

public class Euler778 {

    static final long MOD = 1000000009L;
    static final int D = 10;

    static long[][] matMul(long[][] A, long[][] B) {
        long[][] C = new long[D][D];
        for (int i = 0; i < D; ++i) {
            for (int k = 0; k < D; ++k) {
                long aik = A[i][k];
                if (aik == 0)
                    continue;
                for (int j = 0; j < D; ++j) {
                    C[i][j] = (C[i][j] + aik * B[k][j]) % MOD;
                }
            }
        }
        return C;
    }

    static long[][] matPow(long[][] base, long exp) {
        long[][] result = new long[D][D];
        for (int i = 0; i < D; ++i) {
            result[i][i] = 1;
        }

        while (exp > 0) {
            if ((exp & 1) == 1) {
                result = matMul(base, result);
            }
            base = matMul(base, base);
            exp >>= 1;
        }
        return result;
    }

    static long[] digitCounts(long M, long p10) {
        long N = M + 1;
        long block = 10L * p10;
        long full = N / block;
        long rem = N % block;

        long[] c = new long[D];
        for (int d = 0; d < D; ++d) {
            long cnt = full * p10;
            long shift = d * p10;
            if (rem > shift) {
                long extra = rem - shift;
                cnt += (extra < p10) ? extra : p10;
            }
            c[d] = cnt % MOD;
        }
        return c;
    }

    static long contributionForPosition(long R, long M, long p10) {
        long[] c = digitCounts(M, p10);

        long[][] T = new long[D][D];
        for (int from = 0; from < D; ++from) {
            for (int d = 0; d < D; ++d) {
                int to = (from * d) % 10;
                T[to][from] = (T[to][from] + c[d]) % MOD;
            }
        }

        long[][] P = matPow(T, R);

        long s = 0;
        for (int r = 0; r < D; ++r) {
            s = (s + r * P[r][1]) % MOD;
        }
        return s;
    }

    static long fMod(long R, long M) {
        long ans = 0;
        long placeMod = 1;

        for (long p10 = 1; p10 <= M; p10 *= 10L) {
            long s = contributionForPosition(R, M, p10);
            ans = (ans + placeMod * s) % MOD;
            placeMod = (placeMod * 10L) % MOD;
        }

        return ans;
    }

    public static String solve() {
        return Long.toString(fMod(234567, 765432));
    }

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