Detailed mathematical approach

Problem Summary

Let \(k\) be the number of decimal digits of \(d\). We want to count

$F(n,d)=\#\left\{x\in \mathbb{Z}_{>0}:x\mid n!,\ x\equiv d \pmod{10^k}\right\},$

with the final result reported modulo \(10^{16}+61\). For Project Euler 474 the concrete input is \(n=10^6\) and \(d=65432\). The difficulty is that \(n!\) has an enormous number of divisors, so the computation must work with prime exponents and modular residues rather than explicit divisor generation.

Mathematical Approach

Every divisor of \(n!\) is determined by choosing one exponent for each prime \(p\le n\). The implementation uses that representation, then separates the forced powers of \(2\) and \(5\) coming from the target suffix, and finally runs a dynamic program on the unit group modulo a reduced power of \(10\).

Step 1: Describe Divisors of \(n!\) by Prime Exponents

If

$n!=\prod_{p\le n} p^{e_p},\qquad e_p=v_p(n!)=\sum_{m\ge 1}\left\lfloor\frac{n}{p^m}\right\rfloor,$

then every divisor \(x\mid n!\) has the form

$x=\prod_{p\le n} p^{a_p},\qquad 0\le a_p\le e_p.$

So the problem is equivalent to counting exponent vectors \((a_p)\) that place the resulting divisor in one prescribed residue class modulo \(10^k\).

Step 2: Remove the Exact Powers of \(2\) and \(5\) Forced by the Suffix

Write

$d=2^{\alpha}5^{\beta}\tau,\qquad \gcd(\tau,10)=1,$

where \(\alpha=v_2(d)\) and \(\beta=v_5(d)\). In the main regime used for the actual Euler input, both \(\alpha\lt k\) and \(\beta\lt k\). Define

$f=2^{\alpha}5^{\beta},\qquad M=\frac{10^k}{f}=2^{k-\alpha}5^{k-\beta},\qquad t=\frac{d}{f}.$

If a divisor \(x\) satisfies \(x\equiv d \pmod{10^k}\), then \(x\) must contain exactly \(\alpha\) factors of \(2\) and exactly \(\beta\) factors of \(5\). Dividing the congruence by \(f\) gives

$\frac{x}{f}\equiv t \pmod{M},$

and \(t\) is coprime to \(10\), so \(x/f\) is automatically a unit modulo \(M\). Therefore the exponents of \(2\) and \(5\) are fixed, and only the primes \(p\ne 2,5\) remain free.

If \(\alpha>v_2(n!)\) or \(\beta>v_5(n!)\), the answer is immediately \(0\). The implementations also keep a separate tiny brute-force branch for the checkpoint \((n,d)=(12,12)\), because there \(\alpha=k\) and the unit reduction above does not apply.

Step 3: Dynamic Programming on the Unit Residues Modulo \(M\)

After factoring out the fixed part \(f\), the remaining unit part of a divisor is

$u=\prod_{\substack{p\le n\\ p\ne 2,5}} p^{a_p}\pmod{M}.$

Only residues coprime to \(M\) can occur, so the state space has size

$U=\varphi(M).$

Define \(A(r)\) as the number of choices of exponents for the primes processed so far that produce residue \(r\) modulo \(M\). Initially only the empty product is present, so

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

When processing one prime \(p\ne 2,5\) with exponent range \(0\le j\le e_p\), the transition is

$A_{\mathrm{new}}(x)=\sum_{j=0}^{e_p} A_{\mathrm{old}}(x p^{-j}),$

where the inverse is taken modulo \(M\). This is correct because choosing exponent \(j\) for \(p\) multiplies the current residue by \(p^j\).

Step 4: Use Cycle Decomposition to Make Each Prime Update Fast

Since \(\gcd(p,M)=1\), multiplication by \(p\) permutes the unit residues modulo \(M\). Hence the unit set splits into disjoint cycles

$u_0,\ u_1,\ \dots,\ u_{L-1},\qquad u_{i+1}\equiv p\,u_i \pmod{M}.$

On one cycle, the transition becomes a cyclic convolution:

$A_{\mathrm{new}}(u_i)=\sum_{j=0}^{e_p} A_{\mathrm{old}}(u_{i-j}),$

with indices interpreted modulo \(L\). Write

$e_p+1=qL+r,\qquad 0\le r\lt L.$

Then each new value contains \(q\) complete wraps around the cycle plus an extra window of length \(r\):

$A_{\mathrm{new}}(u_i)=q\sum_{m=0}^{L-1}A_{\mathrm{old}}(u_m)+\sum_{j=0}^{r-1}A_{\mathrm{old}}(u_{i-j}).$

The first term is constant on the whole cycle, and the second can be updated with a sliding window. That reduces the work for one prime from \(O(e_pL)\) on a cycle to \(O(L)\).

Step 5: Read the Target Residue

After all primes \(p\ne 2,5\) have been processed, the desired count is simply the state at the target unit residue \(t\):

$F(n,d)\equiv A(t)\pmod{10^{16}+61}.$

All arithmetic in the dynamic program is performed modulo \(10^{16}+61\), while the residue classes themselves live modulo \(M\).

Worked Example: The Actual Suffix \(d=65432\)

Here \(k=5\) and

$65432=2^3\cdot 8179.$

So \(\alpha=3\), \(\beta=0\), \(f=8\), and

$M=\frac{10^5}{8}=12500,\qquad t=\frac{65432}{8}=8179.$

Because \(\gcd(8179,12500)=1\), the main unit-case algorithm applies directly. The counted divisors are exactly those divisors of \(n!\) whose exponent of \(2\) is \(3\), whose exponent of \(5\) is \(0\), and whose remaining unit part is congruent to \(8179\) modulo \(12500\). The dynamic program therefore runs on

$U=\varphi(12500)=12500\left(1-\frac12\right)\left(1-\frac15\right)=5000$

unit residues instead of on the enormous set of all divisors of \(10^6!\).

How the Code Works

The C++, Python, and Java implementations follow the same numerical strategy. They first sieve all primes up to \(n\), compute \(v_p(n!)\) by Legendre's formula, and determine whether the target suffix falls into the generic unit-case or into the one small checkpoint handled separately by direct divisor enumeration.

In the unit-case, the implementation builds the list of residues modulo \(M\) that are coprime to \(M\), maps each such residue to an array position, and initializes a one-dimensional dynamic-programming table with count \(1\) at residue \(1\).

For each prime \(p\ne 2,5\), it reuses the cycle decomposition of the permutation \(r\mapsto pr \pmod{M}\). Inside each cycle it computes the full-cycle contribution once, then advances the remaining partial sum with a sliding window, so every state in that cycle is updated in constant amortized time.

Because the answer modulus is slightly above \(10^{16}\), the implementations also use overflow-safe modular arithmetic for the running totals. After all odd primes have been processed, the entry corresponding to \(t=d/f\) is the required count modulo \(10^{16}+61\).

Complexity Analysis

Let \(U=\varphi(M)\), and let \(C\) be the number of distinct residue classes \(p\bmod M\) that occur among the primes \(p\le n\) with \(p\ne 2,5\). The prime sieve costs \(O(n\log\log n)\) time and \(O(n)\) memory. Building the cached cycle layouts costs \(O(CU)\) time and \(O(CU)\) memory in the current implementation. Once those layouts are available, each prime transition touches each unit residue only once, so the dynamic-programming work is \(O(\pi(n)\,U)\).

For the actual Euler input, \(M=12500\) and \(U=5000\), which is why this approach is practical even though \(10^6!\) itself has an astronomically large divisor set.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=474
2. Legendre's formula: Wikipedia — Legendre's formula
3. Euler's totient function: Wikipedia — Euler's totient function
4. Multiplicative group of integers modulo \(n\): Wikipedia — Multiplicative group of integers modulo \(n\)
5. Permutations and cycle decomposition: Wikipedia — Permutation