Project Euler Problem 429

Solution

Answer: 98792821

Let

$$n! = \prod_{p \le n} p^{e_p}$$

where

$$e_p = \sum_{k\ge1}\left\lfloor \frac{n}{p^k}\right\rfloor$$

by Legendre’s formula.

A unitary divisor $d$ of $n!$ must satisfy

$$\gcd(d,n!/d)=1.$$

Since the prime factorization of $n!$ is unique, for each prime $p^{e_p}$, a unitary divisor either takes the entire factor $p^{e_p}$ or none of it. Taking only part of the exponent would make the gcd nontrivial.

Thus every unitary divisor has the form

$$d=\prod_{p\le n} p^{\varepsilon_p e_p}, \qquad \varepsilon_p\in{0,1}.$$

We want the sum of squares:

$$S(n!) = \sum_{d \text{ unitary}} d^2.$$

Because the choices are independent across primes, this factors multiplicatively:

$$S(n!) = \prod_{p\le n} \left(1+p^{2e_p}\right).$$

So for $n=100,000,000$,

$$S(100000000!) = \prod_{p\le 100000000} \left(1+p^{2e_p}\right) \pmod{1000000009}.$$

We compute:

Generate all primes $p\le 10^8$ via a sieve.
For each prime, compute $e_p$ using Legendre’s formula.
Multiply

$$(1+p^{2e_p}) \bmod 1000000009$$

into the running product.

Python implementation:

from math import isqrt

MOD = 1_000_000_009
N = 100_000_000

# Sieve of Eratosthenes
is_prime = bytearray(b'\x01') * (N + 1)
is_prime[0:2] = b'\x00\x00'

for i in range(2, isqrt(N) + 1):
    if is_prime[i]:
        is_prime[i * i : N + 1 : i] = (
            b'\x00' * (((N - i * i) // i) + 1)
        )

primes = [p for p in range(2, N + 1) if is_prime[p]]

result = 1

for p in primes:
    # Legendre exponent of p in N!
    e = 0
    m = N
    while m:
        m //= p
        e += m

    result = (result * (1 + pow(p, 2 * e, MOD))) % MOD

print(result)

Small check:

For $4!=24=2^3\cdot 3$,

$$S(4!) = (1+2^{6})(1+3^2) = (1+64)(1+9) = 65\cdot 10 = 650,$$

matching the example.

Running the full computation for $100000000!$ gives:

Answer: 98792821