Problem 627: Counting Products

Project Euler Problem 627 Solution

EulerSolve provides an optimized solution for Project Euler Problem 627, Counting Products, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For fixed integers $m$ and $n$, let $F(m,n)$ denote the number of distinct integers that can be written in the form $$a_1a_2\cdots a_n,\qquad 1\le a_i\le m.$$ The target case is $F(30,10001)$ modulo $10^9+7$. Directly exploring all $30^{10001}$ sequences is impossible, and even storing all distinct products becomes unmanageable very quickly. The solution therefore uses a closed generating-function identity specialized to $m=30$. Mathematical Approach The key observation is that equal products have equal prime-exponent data. Once the problem is translated into exponent vectors, the sequence $F(30,n)$ is encoded by a rational generating function, and the coefficient of $x^n$ can be extracted explicitly. Step 1: Replace Products by Prime-Exponent Vectors Every integer $a$ with $1\le a\le 30$ has a unique factorization over the ten primes not exceeding $30$: $$2,3,5,7,11,13,17,19,23,29.$$ Write $$a=\prod_{p\le 30} p^{e_p(a)}.$$ Then the product $a_1a_2\cdots a_n$ is determined by the sum of these exponent vectors: $$E(a_1,\dots,a_n)=\sum_{j=1}^{n}\bigl(e_p(a_j)\bigr)_{p\le 30}.$$ Two products are equal if and only if their exponent sums are equal. Therefore $F(30,n)$ is exactly the number of lattice points in $\mathbb{Z}_{\ge 0}^{10}$ that can be reached as sums of $n$ generator vectors coming from the numbers $1,2,\dots,30$....

Detailed mathematical approach

Problem Summary

For fixed integers $m$ and $n$, let $F(m,n)$ denote the number of distinct integers that can be written in the form

$$a_1a_2\cdots a_n,\qquad 1\le a_i\le m.$$

The target case is $F(30,10001)$ modulo $10^9+7$. Directly exploring all $30^{10001}$ sequences is impossible, and even storing all distinct products becomes unmanageable very quickly. The solution therefore uses a closed generating-function identity specialized to $m=30$.

Mathematical Approach

The key observation is that equal products have equal prime-exponent data. Once the problem is translated into exponent vectors, the sequence $F(30,n)$ is encoded by a rational generating function, and the coefficient of $x^n$ can be extracted explicitly.

Step 1: Replace Products by Prime-Exponent Vectors

Every integer $a$ with $1\le a\le 30$ has a unique factorization over the ten primes not exceeding $30$:

$$2,3,5,7,11,13,17,19,23,29.$$

Write

$$a=\prod_{p\le 30} p^{e_p(a)}.$$

Then the product $a_1a_2\cdots a_n$ is determined by the sum of these exponent vectors:

$$E(a_1,\dots,a_n)=\sum_{j=1}^{n}\bigl(e_p(a_j)\bigr)_{p\le 30}.$$

Two products are equal if and only if their exponent sums are equal. Therefore $F(30,n)$ is exactly the number of lattice points in $\mathbb{Z}_{\ge 0}^{10}$ that can be reached as sums of $n$ generator vectors coming from the numbers $1,2,\dots,30$. Multiplication has been converted into addition.

Step 2: Package All Lengths into One Generating Function

For the fixed bound $m=30$, those reachable exponent vectors form a finitely generated commutative semigroup. Counting distinct products of length $n$ is therefore equivalent to taking the coefficient of $x^n$ in the corresponding Hilbert series.

The central structural identity used by the implementation is

$$\sum_{n\ge 0} F(30,n)x^n=\frac{1+19x+33x^2+6x^3}{(1-x)^{11}}.$$

The denominator $(1-x)^{11}$ shows that the sequence behaves like a polynomial sequence of degree $10$, while the short numerator provides the exact initial corrections needed for the fixed generator set $\{1,2,\dots,30\}$.

Step 3: Extract Coefficients with the Binomial Expansion

Use the standard series identity

$$\frac{1}{(1-x)^{11}}=\sum_{t\ge 0}\binom{t+10}{10}x^t.$$

Multiplying by the numerator shifts the indices of those coefficients, so for every $n\ge 0$ we get

$$F(30,n)=\binom{n+10}{10}+19\binom{n+9}{10}+33\binom{n+8}{10}+6\binom{n+7}{10}.$$

If the upper index of a binomial coefficient is smaller than $10$, that term is interpreted as $0$. This makes the formula valid uniformly, including the small cases $n=0,1,2$.

Step 4: Why This Solves the Target Instance

Once the generating function is known, the difficult combinatorics is over. For $n=10001$, the desired count is just the weighted sum

$$F(30,10001)=\binom{10011}{10}+19\binom{10010}{10}+33\binom{10009}{10}+6\binom{10008}{10},$$

followed by reduction modulo $10^9+7$. Instead of tracking a gigantic set of products, we only need four combinatorial terms.

Step 5: Worked Example

For $n=2$, the formula gives

$$F(30,2)=\binom{12}{10}+19\binom{11}{10}+33\binom{10}{10}+6\binom{9}{10}.$$

Since $\binom{9}{10}=0$, this becomes

$$F(30,2)=66+19\cdot 11+33\cdot 1=66+209+33=308.$$

For $n=3$, we similarly obtain

$$F(30,3)=\binom{13}{10}+19\binom{12}{10}+33\binom{11}{10}+6\binom{10}{10}=286+1254+363+6=1909.$$

These values agree with the explicit small-case checks used by the implementation, so they serve as practical confirmation that the closed form has been applied correctly.

How the Code Works

The C++, Python, and Java implementations all evaluate the same closed formula. They precompute factorials and inverse factorials up to $n+10$ modulo $10^9+7$, which allows each binomial coefficient $\binom{r}{10}$ to be recovered in constant time.

Because $10^9+7$ is prime, the inverse factorial table is obtained from one modular inverse and a backward sweep, using Fermat's little theorem:

$$a^{-1}\equiv a^{10^9+5}\pmod{10^9+7}.$$

After that preprocessing, the implementation forms the weighted sum of the four shifted binomial terms with coefficients $1$, $19$, $33$, and $6$. The C++ implementation also verifies several small exact cases by brute force, including $F(9,2)=36$, $F(30,2)=308$, $F(30,3)=1909$, and $F(30,4)=8679$, before evaluating the large target input.

Complexity Analysis

Let the requested length be $n$. Building factorial and inverse-factorial tables up to $n+10$ costs $O(n)$ time and $O(n)$ memory. The final evaluation uses only four binomial coefficients and a few modular multiplications, so it is $O(1)$ after preprocessing. This is exponentially better than enumerating all $30^n$ sequences or storing all distinct products directly.

Footnotes and References

Problem page: https://projecteuler.net/problem=627
Generating functions: Wikipedia — Generating function
Hilbert series: Wikipedia — Hilbert series
Fundamental theorem of arithmetic: Wikipedia — Fundamental theorem of arithmetic
Binomial coefficient: Wikipedia — Binomial coefficient

Problem 627 source code

C++

#include <cassert>
#include <cstdint>
#include <iostream>
#include <unordered_set>
#include <vector>

using i64 = long long;
using u64 = unsigned long long;

static constexpr int MOD = 1'000'000'007;

static inline int mod_add(int a, int b) {
    int s = a + b;
    if (s >= MOD) s -= MOD;
    if (s < 0) s += MOD;
    return s;
}

static inline int mod_mul(i64 a, i64 b) { return (int)((a % MOD) * (b % MOD) % MOD); }

static int mod_pow(int a, int e) {
    long long r = 1, x = a;
    while (e > 0) {
        if (e & 1) r = (r * x) % MOD;
        x = (x * x) % MOD;
        e >>= 1;
    }
    return (int)r;
}

struct Comb {
    std::vector<int> fact, ifact;

    explicit Comb(int n) : fact(n + 1, 1), ifact(n + 1, 1) {
        for (int i = 1; i <= n; ++i) fact[i] = mod_mul(fact[i - 1], i);
        ifact[n] = mod_pow(fact[n], MOD - 2);
        for (int i = n; i >= 1; --i) ifact[i - 1] = mod_mul(ifact[i], i);
    }

    int C(int n, int k) const {
        if (k < 0 || k > n) return 0;
        return mod_mul(fact[n], mod_mul(ifact[k], ifact[n - k]));
    }
};

static int F30(int n, const Comb &cb) {
    int ans = 0;
    ans = mod_add(ans, cb.C(n + 10, 10));
    if (n >= 1) ans = mod_add(ans, mod_mul(19, cb.C(n + 9, 10)));
    if (n >= 2) ans = mod_add(ans, mod_mul(33, cb.C(n + 8, 10)));
    if (n >= 3) ans = mod_add(ans, mod_mul(6, cb.C(n + 7, 10)));
    return ans;
}

static int F_bruteforce(int m, int n) {
    std::unordered_set<u64> seen;
    seen.reserve(1u << 20);
    std::vector<u64> cur{1};
    for (int step = 0; step < n; ++step) {
        std::unordered_set<u64> nxt;
        nxt.reserve(cur.size() * (size_t)m);
        for (u64 x : cur) {
            for (int a = 1; a <= m; ++a) nxt.insert(x * (u64)a);
        }
        cur.assign(nxt.begin(), nxt.end());
    }
    return (int)cur.size();
}

int main() {
    const int N = 10'001;
    Comb cb(N + 10);

    assert(F_bruteforce(9, 2) == 36);
    assert(F_bruteforce(30, 2) == 308);
    assert(F_bruteforce(30, 3) == 1909);
    assert(F_bruteforce(30, 4) == 8679);

    assert(F30(0, cb) == 1);
    assert(F30(1, cb) == 30);
    assert(F30(2, cb) == 308);
    assert(F30(3, cb) == 1909);
    assert(F30(4, cb) == 8679);
    assert(F30(11, cb) == 7'174'102);
    assert(F30(15, cb) == 82'166'678);

    std::cout << F30(N, cb) << "\n";
    return 0;
}

Python

MOD = 1000000007

def mod_pow(a, e):
    return pow(a, e, MOD)

class Comb:
    def __init__(self, n):
        self.fact = [1] * (n + 1)
        self.ifact = [1] * (n + 1)
        for i in range(1, n + 1):
            self.fact[i] = (self.fact[i - 1] * i) % MOD
        self.ifact[n] = mod_pow(self.fact[n], MOD - 2)
        for i in range(n, 0, -1):
            self.ifact[i - 1] = (self.ifact[i] * i) % MOD
            
    def C(self, n, k):
        if k < 0 or k > n: return 0
        return (self.fact[n] * self.ifact[k] % MOD) * self.ifact[n - k] % MOD

def F30(n, cb):
    ans = 0
    ans = (ans + cb.C(n + 10, 10)) % MOD
    if n >= 1:
        ans = (ans + 19 * cb.C(n + 9, 10)) % MOD
    if n >= 2:
        ans = (ans + 33 * cb.C(n + 8, 10)) % MOD
    if n >= 3:
        ans = (ans + 6 * cb.C(n + 7, 10)) % MOD
    return ans

def solve():
    N = 10001
    cb = Comb(N + 10)
    return str(F30(N, cb))

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

Java

public class Euler627 {
    static final int MOD = 1000000007;

    static int modAdd(int a, int b) {
        int s = a + b;
        if (s >= MOD)
            s -= MOD;
        if (s < 0)
            s += MOD;
        return s;
    }

    static int modMul(long a, long b) {
        return (int) ((a % MOD) * (b % MOD) % MOD);
    }

    static int modPow(int a, int e) {
        long r = 1, x = a;
        while (e > 0) {
            if ((e & 1) == 1)
                r = (r * x) % MOD;
            x = (x * x) % MOD;
            e >>= 1;
        }
        return (int) r;
    }

    static class Comb {
        int[] fact;
        int[] ifact;

        Comb(int n) {
            fact = new int[n + 1];
            ifact = new int[n + 1];
            fact[0] = 1;
            ifact[0] = 1;

            for (int i = 1; i <= n; i++)
                fact[i] = modMul(fact[i - 1], i);
            ifact[n] = modPow(fact[n], MOD - 2);
            for (int i = n; i >= 1; i--)
                ifact[i - 1] = modMul(ifact[i], i);
        }

        int C(int n, int k) {
            if (k < 0 || k > n)
                return 0;
            return modMul(fact[n], modMul(ifact[k], ifact[n - k]));
        }
    }

    static int F30(int n, Comb cb) {
        int ans = 0;
        ans = modAdd(ans, cb.C(n + 10, 10));
        if (n >= 1)
            ans = modAdd(ans, modMul(19, cb.C(n + 9, 10)));
        if (n >= 2)
            ans = modAdd(ans, modMul(33, cb.C(n + 8, 10)));
        if (n >= 3)
            ans = modAdd(ans, modMul(6, cb.C(n + 7, 10)));
        return ans;
    }

    public static String solve() {
        int N = 10001;
        Comb cb = new Comb(N + 10);
        return Integer.toString(F30(N, cb));
    }

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