Detailed mathematical approach

Problem Summary

Let \(W(n,k)\) denote the number of ways to write \(n\) as a product of \(k\) distinct positive integers, with order ignored. The target instance in this problem is

$W(10000!,30)\pmod{10^9+7}.$

A direct search over divisors or over candidate factor sets is hopelessly large. The implementation therefore works with the prime-exponent vector of \(n\), converts the distinct-factor condition into a generating-function coefficient, and evaluates that coefficient through a sum over integer partitions of \(k\).

Mathematical Approach

Write the prime factorization of \(n\) as

$n=\prod_p p^{e_p}.$

For the factorial input \(n=N!\), the exponents are obtained from Legendre's formula:

$e_p=\sum_{j\ge 1}\left\lfloor\frac{N}{p^j}\right\rfloor.$

Any factor in a valid decomposition can be written in the form

$a_i=\prod_p p^{\alpha_{i,p}},\qquad 0\le \alpha_{i,p}\le e_p,$

and the condition \(a_1a_2\cdots a_k=n\) is equivalent to

$\alpha_{1,p}+\alpha_{2,p}+\cdots+\alpha_{k,p}=e_p \qquad \text{for every prime } p.$

Step 1: Encode Unordered Distinct Factor Selections

For each divisor \(d\mid n\), introduce the monomial

$m(d)=\prod_p x_p^{v_p(d)}.$

Now consider the product

$F(t,\mathbf{x})=\prod_{d\mid n}\left(1+t\,m(d)\right).$

Selecting the term \(t\,m(d)\) means that the divisor \(d\) is chosen; selecting \(1\) means that it is skipped. Because each divisor appears only once in the product, it can be selected at most once, so the distinctness condition is built in automatically.

The coefficient

$[t^k\prod_p x_p^{e_p}]\,F(t,\mathbf{x})$

therefore counts unordered sets of \(k\) distinct divisors whose product is exactly \(n\). Hence

$W(n,k)=[t^k\prod_p x_p^{e_p}]\,F(t,\mathbf{x}).$

Step 2: Replace the Subset Product by Symmetric Polynomials

If we denote by \(\mathcal{M}\) the set of all divisor monomials \(m(d)\), then

$F(t,\mathbf{x})=\sum_{r\ge 0} e_r(\mathcal{M})\,t^r,$

where \(e_r\) is the \(r\)-th elementary symmetric polynomial. So the problem is to extract the coefficient of \(\prod_p x_p^{e_p}\) inside \(e_k(\mathcal{M})\).

A standard identity expands \(e_k\) in terms of power sums:

$e_k=\sum_{\lambda\vdash k}\frac{(-1)^{k-\ell(\lambda)}}{z_\lambda}\,p_\lambda,$

where \(\lambda=(\lambda_1,\dots,\lambda_r)\) is an integer partition of \(k\), \(\ell(\lambda)=r\), \(p_\lambda=\prod_{i=1}^r p_{\lambda_i}\), and

$z_\lambda=\prod_v v^{m_v}m_v!.$

Here \(m_v\) denotes how many times the part value \(v\) occurs in \(\lambda\). This identity is the reason the implementation iterates over partitions of \(k\) instead of over subsets of divisors.

Step 3: The Partition Coefficient Used in the Program

The coefficient attached to \(\lambda\) can be rewritten as

$\frac{(-1)^{k-\ell(\lambda)}}{z_\lambda} =\left(\prod_{i=1}^r\frac{(-1)^{\lambda_i-1}}{\lambda_i}\right)\left(\prod_v\frac{1}{m_v!}\right).$

So every partition contributes a rational weight determined by the sizes of its parts and by the multiplicities of repeated parts. In the final modular computation, those divisions are implemented with modular inverses modulo

$M=10^9+7.$

This matches the coefficient pattern used by the C++, Python, and Java implementations.

Step 4: Evaluate One Partition Prime by Prime

For a positive integer \(j\), the power sum is

$p_j=\sum_{d\mid n} m(d)^j.$

Since divisors are described independently prime by prime, this factorizes as

$p_j=\prod_p \sum_{a=0}^{e_p} x_p^{ja}.$

Now fix a partition \(\lambda=(\lambda_1,\dots,\lambda_r)\). For a single prime exponent \(e\), the needed one-variable coefficient is

$d_\lambda(e)=[x^e]\prod_{i=1}^r\frac{1}{1-x^{\lambda_i}}.$

This is safe because only terms of degree at most \(e\) can contribute to \([x^e]\). Equivalently, \(d_\lambda(e)\) is the number of nonnegative integer solutions of

$\lambda_1 b_1+\lambda_2 b_2+\cdots+\lambda_r b_r=e.$

The implementation computes these coefficients with an unbounded-knapsack dynamic program up to the largest exponent that occurs in the factorization of \(n\).

Step 5: Group Equal Prime Exponents

Many primes share the same exponent in \(n\). Let

$f_e=\#\{p:e_p=e\}.$

For a fixed partition \(\lambda\), all primes with the same exponent contribute the same coefficient, so the total partition weight becomes

$S(\lambda)=\prod_e d_\lambda(e)^{f_e}.$

Therefore the full answer is

$W(n,k)=\sum_{\lambda\vdash k}\left(\frac{(-1)^{k-\ell(\lambda)}}{z_\lambda}\right)S(\lambda)\pmod{M}.$

For the target \(n=N!\), the exponent profile \(e_p\) comes from Legendre's formula, and grouping equal values avoids repeating the same DP lookup for many different primes.

Step 6: Worked Example \(W(144,4)=7\)

Take

$144=2^4\cdot 3^2.$

So the exponent multiset is \(\{4,2\}\), which means \(f_2=1\) and \(f_4=1\). The partitions of \(4\) are

$1+1+1+1,\qquad 2+1+1,\qquad 2+2,\qquad 3+1,\qquad 4.$

For each partition we evaluate its coefficient and the two needed values \(d_\lambda(2)\) and \(d_\lambda(4)\):

$\begin{aligned} \lambda=1+1+1+1&:\quad \text{coeff}=\frac{1}{24},\quad d_\lambda(2)=10,\quad d_\lambda(4)=35,\quad \text{term}=\frac{350}{24},\\ \lambda=2+1+1&:\quad \text{coeff}=-\frac{1}{4},\quad d_\lambda(2)=4,\quad d_\lambda(4)=9,\quad \text{term}=-9,\\ \lambda=2+2&:\quad \text{coeff}=\frac{1}{8},\quad d_\lambda(2)=2,\quad d_\lambda(4)=3,\quad \text{term}=\frac{3}{4},\\ \lambda=3+1&:\quad \text{coeff}=\frac{1}{3},\quad d_\lambda(2)=1,\quad d_\lambda(4)=2,\quad \text{term}=\frac{2}{3},\\ \lambda=4&:\quad \text{coeff}=-\frac{1}{4},\quad d_\lambda(2)=0,\quad d_\lambda(4)=1,\quad \text{term}=0. \end{aligned}$

Summing the five terms gives

$\frac{350}{24}-9+\frac{3}{4}+\frac{2}{3}=7.$

So \(W(144,4)=7\), exactly matching the checkpoint used by the implementation.

How the Code Works

For the factorial target, the implementation first generates the relevant primes and computes each exponent \(e_p\) by Legendre's formula. It then compresses that data into the frequency table \(f_e\), because only the exponent values and how often they occur matter in the partition formula.

Next, it generates all integer partitions of \(k\). For each partition, it computes the modular version of the coefficient \(\frac{(-1)^{k-\ell(\lambda)}}{z_\lambda}\) using fast exponentiation for modular inverses and precomputed inverse factorials for repeated part multiplicities.

After that, the implementation runs an unbounded-knapsack DP up to the maximum exponent \(E=\max e_p\). The DP value at index \(e\) is precisely \(d_\lambda(e)\), the coefficient of \(x^e\) in \(\prod_i(1-x^{\lambda_i})^{-1}\). Using the frequency table \(f_e\), it multiplies the relevant DP values, raises them to the needed powers, and adds the weighted contribution of the partition to the total.

Because different partitions are independent, the C++, Python, and Java implementations process disjoint partition ranges in parallel and combine the partial sums modulo \(10^9+7\).

Complexity Analysis

Let \(p(k)\) be the number of integer partitions of \(k\), and let \(E=\max_p e_p\). For the factorial input \(n=N!\), generating primes and evaluating Legendre's formula with the straightforward sieve-based approach costs \(O(N\log\log N)\) time and \(O(N)\) memory.

For a fixed partition \(\lambda\), the dynamic program runs to degree \(E\) and performs one complete-knapsack pass per part, so its cost is \(O(\ell(\lambda)E)\subseteq O(kE)\). Summing over all partitions gives total combinatorial work \(O(p(k)kE)\). The working memory is \(O(E)\) per worker, plus storage for the partition list and the exponent-frequency table. Parallelism reduces wall-clock time but not the asymptotic arithmetic count.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=495
2. Legendre's formula: Wikipedia - Legendre's formula
3. Integer partition: Wikipedia - Integer partition
4. Elementary symmetric polynomial: Wikipedia - Elementary symmetric polynomial
5. Generating function: Wikipedia - Generating function