Detailed mathematical approach

Problem Summary

The quantity \(F(n)\) counts restricted factorisations of \(n!\) of the shape

$n!=b_1\,b_2^2\,b_3^2\,b_4^3\,b_5^3\,b_6^3\,b_7^4\,b_8^4\,b_9^4\,b_{10}^4.$

The computation is organized in two stages. First, the ten slots are treated as labeled so that equal-base collisions can be handled systematically by inclusion-exclusion. Second, the labels inside the equal-exponent groups are forgotten, because swapping the two square slots, the three cube slots, or the four fourth-power slots does not change the underlying factorisation. The goal is to evaluate \(F(10^6)\bmod(10^9+7)\).

Mathematical Approach

Let the exponent pattern be

$\alpha=(\alpha_1,\dots,\alpha_{10})=(1,2,2,3,3,3,4,4,4,4).$

If we temporarily label the bases as \(b_1,\dots,b_{10}\), then the only real difficulty is that some of these bases may coincide. The implementations resolve that by working on the partition lattice of the ten labels and by treating each prime of \(n!\) independently.

Step 1: Compute the prime exponents of \(n!\)

For each prime \(p\le n\), Legendre's formula gives

$e_p=v_p(n!)=\sum_{k\ge 1}\left\lfloor\frac{n}{p^k}\right\rfloor.$

Prime powers are independent, so the global count factors over the multiset of values \(\{e_p\}\). It is therefore enough to store the histogram

$c_e=\#\{p\le n:\ v_p(n!)=e\}.$

Once the numbers \(c_e\) are known, every partition profile contributes through repeated powers of the same local coefficient \(g(e)\).

Step 2: Encode equal-base collisions by set partitions

Suppose several labeled slots share the same base. Their collision pattern is a set partition \(\pi\) of \(\{1,\dots,10\}\). Every block \(B\in\pi\) represents one actual base, and the exponents attached to that base add up to

$S_B=\sum_{i\in B}\alpha_i.$

So after choosing \(\pi\), the original labels only matter through the multiset of block sums \(\{S_B\}\). Two different partitions that produce the same sorted block-sum vector have the same local generating function, so they can be merged by adding their Möbius weights. The implementations do exactly that.

Step 3: Count one prime with a generating function

Fix a prime \(p\) with exponent \(e=e_p\). If block \(B\) receives \(u_B\ge 0\) copies of \(p\) in its common base, then the exponent balance is

$\sum_{B\in\pi} S_Bu_B=e.$

The number of nonnegative solutions is the coefficient of

$G_\pi(x)=\prod_{B\in\pi}\frac{1}{1-x^{S_B}}=\sum_{m\ge 0} g_\pi(m)x^m.$

Therefore \(g_\pi(e)\) is exactly the number of ways a single prime of valuation \(e\) can be distributed across the bases while respecting the collision pattern \(\pi\).

Step 4: Recover distinct labeled bases by Möbius inversion

Counting every partition \(\pi\) would allow equal bases. To force the ten labeled bases to be distinct, we use Möbius inversion on the partition lattice. Its weight is

$\mu(\pi)=(-1)^{10-|\pi|}\prod_{B\in\pi}(|B|-1)!.$

Because primes are independent, the total labeled contribution of \(\pi\) is

$W_\pi(n)=\prod_{p\le n} g_\pi(e_p)=\prod_{e\ge 1} g_\pi(e)^{c_e}.$

Hence the number of factorizations with fully labeled and pairwise distinct bases is

$L(n)=\sum_{\pi}\mu(\pi)\,W_\pi(n).$

Step 5: Forget relabelings inside equal-exponent groups

The labeled model distinguishes the two exponent-\(2\) slots, the three exponent-\(3\) slots, and the four exponent-\(4\) slots. The original problem does not. Every genuine restricted factorisation is therefore counted

$2!\,3!\,4!=288$

times inside \(L(n)\). The desired answer is

$F(n)=\frac{L(n)}{288}\pmod{10^9+7}.$

Since the modulus is prime, the division is implemented by multiplying with the modular inverse of \(288\).

Step 6: Evaluate \(g_\pi(e)\) quickly for all relevant exponents

For small exponents, \(g_\pi(e)\) is computed directly by unbounded knapsack on the block sums \(S_B\). This is exactly coefficient extraction for \(G_\pi(x)\).

For large exponents, the denominator of \(G_\pi(x)\) has degree at most

$1+2+2+3+3+3+4+4+4+4=30,$

so if

$Q_\pi(x)=\prod_{B\in\pi}(1-x^{S_B})=1+q_1x+\cdots+q_{30}x^{30},$

then the coefficients satisfy the linear recurrence

$g_\pi(m)=-\sum_{j=1}^{30}q_j\,g_\pi(m-j)\qquad(m\ge 30).$

This fixed order is the key reason the implementation can jump to distant coefficients quickly instead of extending the dynamic program all the way to \(v_2(n!)\).

Worked Example: the partition with no collisions

If every labeled base is distinct, all ten blocks are singletons and the local series becomes

$G(x)=\frac{1}{(1-x)(1-x^2)^2(1-x^3)^3(1-x^4)^4}.$

Now read off a few coefficients. We have \(g(1)=1\), because an exponent \(1\) can only come from the first slot. We also have \(g(2)=3\), because exponent \(2\) can be formed either as \(2\cdot 1\) in the first slot or as one copy of either of the two square slots.

This tiny example illustrates two implementation ideas at once. First, the local coefficients are genuine coin-change counts for the block sums. Second, if a partition profile had all block sums divisible by some \(d>1\), then \(g_\pi(1)=0\). For the target input \(n=10^6\), primes with valuation \(1\) do occur in \(n!\), so such profiles can be discarded immediately.

How the Code Works

The C++, Python, and Java implementations follow the same mathematical pipeline. They first enumerate set partitions of the ten exponent slots, convert each partition to its sorted block-sum vector, merge profiles with identical vectors, and store the corresponding Möbius weight. They also keep the gcd of each profile so that obviously impossible profiles can be pruned when exponent \(1\) appears in the histogram.

Next, the implementation sieves primes up to \(n\), applies Legendre's formula to every prime, and compresses the result into the frequency table \(c_e\). For each surviving profile it computes all small coefficients \(g_\pi(e)\) by dynamic programming, then uses the degree-\(30\) recurrence to evaluate larger exponents. The contribution \(\mu(\pi)\prod_e g_\pi(e)^{c_e}\) is accumulated modulo \(10^9+7\), and the final sum is multiplied by the modular inverse of \(288\). The C++ version parallelizes independent profile evaluations, the Java version keeps the same arithmetic sequentially, and the Python version delegates to the same C++ core.

Complexity Analysis

The partition-lattice work depends only on the fixed exponent multiset \((1,2,2,3,3,3,4,4,4,4)\), so it is constant-size preprocessing for this problem. After that, building the prime table up to \(n\) costs \(O(n\log\log n)\) time and \(O(n)\) memory. Let \(P\) be the number of distinct block-sum profiles after aggregation; \(P\) is again a problem-specific constant. The dynamic-programming phase for each profile only goes up to about \(\sqrt n\), while larger coefficients are recovered from an order-\(30\) recurrence in logarithmic time. Consequently the overall asymptotic cost is dominated by the sieve, with \(O(n\log\log n)\) time and \(O(n)\) memory, and the partition-profile work contributes a manageable constant factor.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=636
2. Legendre's formula: Wikipedia - Legendre's formula
3. Set partitions and partition lattices: Wikipedia - Partition of a set
4. Inclusion-exclusion principle: Wikipedia - Inclusion-exclusion principle
5. Generating functions: Wikipedia - Generating function
6. Linear recurrences with constant coefficients: Wikipedia - Linear recurrence with constant coefficients