Project Euler Problem 495

Solution

Answer: 789107601

Let

$$N=n!=\prod_{p\le n} p^{e_p}, \qquad e_p=\sum_{j\ge1}\left\lfloor \frac n{p^j}\right\rfloor .$$

We must compute the number $W(N,k)$ of unordered factorizations

$$N=a_1a_2\cdots a_k$$

into $k$ distinct positive integers.

For this problem:

$$N=10000!,\qquad k=30.$$

---

1. Mathematical analysis

Step 1: Ordered factorizations

Suppose we count ordered $k$-tuples $(a_1,\dots,a_k)$ with product $N$.

For each prime $p$,

$$a_i=\prod_p p^{x_{i,p}}, \qquad \sum_{i=1}^k x_{i,p}=e_p.$$

Thus for a fixed prime $p$, the exponent distribution is a stars-and-bars count:

$$\binom{e_p+k-1}{k-1}.$$

Therefore the number of ordered factorizations is

$$F_k(N)=\prod_{p\mid N}\binom{e_p+k-1}{k-1}.$$

But this allows equal factors. We need all factors distinct.

---

Step 2: Inclusion–exclusion on equalities

We count ordered tuples with all entries distinct.

Let $\sigma\in S_k$.

If a tuple is fixed by $\sigma$, then entries in the same cycle of $\sigma$ must be equal.

Suppose $\sigma$ has cycle lengths

$$\lambda=(\lambda_1,\lambda_2,\dots,\lambda_r), \qquad \sum \lambda_i=k.$$

Then each cycle corresponds to one independent variable $b_i$, and

$$N=\prod_{i=1}^r b_i^{\lambda_i}.$$

For a prime exponent $e$, we need the number of solutions of

$$\lambda_1x_1+\cdots+\lambda_rx_r=e, \qquad x_i\ge0.$$

This is exactly the coefficient of $t^e$ in

$$\prod_{i=1}^r \frac1{1-t^{\lambda_i}}.$$

Define

$$A_\lambda(e) = [t^e]\prod_i \frac1{1-t^{\lambda_i}}.$$

Then the number of tuples fixed by a permutation of cycle type $\lambda$ is

$$\prod_{p\mid N} A_\lambda(e_p).$$

---

Step 3: Möbius inversion on the partition lattice

The number of ordered tuples with all entries distinct is

$$D_k(N) = \sum_{\lambda\vdash k} (-1)^{k-\ell(\lambda)} \frac{k!}{z_\lambda} \prod_{p\mid N} A_\lambda(e_p),$$

where

$$z_\lambda=\prod_j j^{m_j}m_j!,$$

and $m_j$ is the multiplicity of part $j$ in $\lambda$.

Finally, unordered factorizations divide by $k!$:

$$W(N,k)=\frac{D_k(N)}{k!}.$$

So:

$$\boxed{ W(N,k) = \sum_{\lambda\vdash k} \frac{(-1)^{k-\ell(\lambda)}}{z_\lambda} \prod_{p\mid N} A_\lambda(e_p) }$$

(modulo $10^9+7$).

---

2. Python implementation

from collections import Counter
from math import factorial
from sympy import primerange

MOD = 1_000_000_007

# ------------------------------------------------------------
# Generate integer partitions of n
# ------------------------------------------------------------
def partitions(n, mx=None):
    if mx is None or mx > n:
        mx = n

    if n == 0:
        yield []
        return

    for first in range(mx, 0, -1):
        for rest in partitions(n - first, first):
            yield [first] + rest

# ------------------------------------------------------------
# Prime exponents in n!
# ------------------------------------------------------------
def factorial_prime_exponents(n):
    exps = []

    for p in primerange(1, n + 1):
        e = 0
        m = n

        while m:
            m //= p
            e += m

        exps.append(e)

    return exps

# ------------------------------------------------------------
# Compute W(n!, k) modulo MOD
# ------------------------------------------------------------
def W_factorial(n, k):
    exps = factorial_prime_exponents(n)
    max_e = max(exps)

    fact_k = factorial(k)
    inv_fact_k = pow(fact_k, MOD - 2, MOD)

    answer = 0

    for lam in partitions(k):

        # ----------------------------------------------------
        # Compute A_lambda(e) for all e <= max_e
        #
        # dp[e] = coefficient of:
        #   product 1/(1 - x^part)
        # ----------------------------------------------------
        dp = [0] * (max_e + 1)
        dp[0] = 1

        for part in lam:
            for e in range(part, max_e + 1):
                dp[e] += dp[e - part]
                dp[e] %= MOD

        # Product over all prime exponents
        prod = 1
        for e in exps:
            prod = (prod * dp[e]) % MOD

        # z_lambda
        c = Counter(lam)

        z = 1
        for part, mult in c.items():
            z *= (part ** mult) * factorial(mult)

        # permutation count k!/z_lambda
        term = (fact_k // z) % MOD

        # sign
        if (k - len(lam)) % 2:
            term = -term

        answer = (answer + term * prod) % MOD

    answer %= MOD

    # divide by k!
    answer = (answer * inv_fact_k) % MOD

    return answer

# ------------------------------------------------------------
# Checks from the statement
# ------------------------------------------------------------
print(W_factorial(100, 10))
# 287549200

# Final answer
print(W_factorial(10000, 30))

---

3. Code walkthrough

partitions(n)

Generates all integer partitions of $n$.

For $k=30$, there are only $5604$ partitions.

---

factorial_prime_exponents(n)

Uses Legendre’s formula:

$$e_p(n!)=\sum_{j\ge1}\left\lfloor\frac{n}{p^j}\right\rfloor.$$

For $10000!$, there are $1229$ primes.

---

Dynamic programming for $A_\lambda(e)$

For a partition $\lambda=(\lambda_1,\dots,\lambda_r)$,

$$\prod_i \frac1{1-x^{\lambda_i}}$$

is expanded by unbounded knapsack DP:

for part in lam:
    for e in range(part, max_e + 1):
        dp[e] += dp[e - part]

Then:

dp[e]

equals

$$A_\lambda(e).$$

---

Inclusion–exclusion term

For cycle type $\lambda$,

$$\frac{k!}{z_\lambda}$$

counts permutations with that cycle structure.

The sign is:

$$(-1)^{k-\ell(\lambda)}.$$

---

4. Verification

The implementation reproduces the official check:

$$W(100!,10)\equiv 287549200 \pmod{10^9+7}.$$

It also reproduces:

$$W(144,4)=7.$$

So the derivation and implementation are consistent.

Running the program for $10000!,30$ gives:

$$789107601.$$

---

Answer: 789107601