Project Euler Problem 956 Solution

Project Euler Problem 956 Statement

Solution

Project Euler Problem 956 overview:

Solution Summary

Definitions:

Let the superfactorial $\mathrm{sf}(n)$ be defined by

$$\mathrm{sf}(n) = \prod_{k=1}^{n} k!$$

Let the super‑duper factorial $\mathrm{sdf}(N)$ be

$$\mathrm{sdf}(N) = \prod_{i=1}^{N} \mathrm{sf}(i) = \prod_{i=1}^{N} \prod_{k=1}^{i} k!$$

The program first computes the prime factorization of $\mathrm{sdf}(N)$:

$$\mathrm{sdf}(N) = \prod_{p \le N} p^{E_p},$$

where $E_p$ is the exponent of the prime $p$ in that product. Then it considers all divisors $d$ of $\mathrm{sdf}(N)$. Each divisor $d$ corresponds to choosing exponents $0 \le e_p \le E_p$ for each prime $p$, so that

$$d = \prod_{p} p^{e_p}$$

It only wants those divisors whose total exponent sum

$$S(d) = \sum_{p} e_p$$

is congruent to 0 mod N.

And instead of just counting such divisors, it wants the sum of their values, modulo $999{,}999{,}001$. So conceptually:

Answer =
sum of all divisors $d$ of sdf(N) such that
$\sum_{p} e_p \equiv 0 \pmod{N}$

Solution Detail

Constants and primes

N, MOD, ROOT = 1000, 999_999_001, 17
A, B = N + 1, N + 2
primes = list(primerange(N + 1))

Here \( N = 1000 \) is the problem’s n.

\( \text{MOD} = 999{,}999{,}001 \) is a prime modulus.

In a prime modulus, the multiplicative group \( (\mathbb{Z}/\text{MOD}\,\mathbb{Z})^{\times} \) is cyclic, so it has primitive roots. \( \text{ROOT} = 17 \) is such a primitive root (pre-chosen).

primes is the list of all primes \( p \le N \), i.e. all primes that can appear in the factorization of \( \mathrm{sdf}(N) \).

\( A, B \) are just \( N + 1, N + 2 \) to make the exponent formula shorter.

Computing the exponents \( E_p \)

exps = []
for p in primes:
    e, t = 0, p
    while t <= N:
        m = N // t
        e += (t*t*m*(m+1)*(2*m+1)//6 - (A+B)*t*m*(m+1)//2 + A*B*m)//2
        t *= p
    exps.append(e)

For each prime \( p \):

• \( t = p, p^2, p^3, \dots \) while \( t \le N \).
• \( m = \left\lfloor \dfrac{N}{t} \right\rfloor \).

For each power \( t = p^j \) we add:

(t*t*m*(m+1)*(2*m+1)//6    # t^2 * S2(m)
 - (A+B)*t*m*(m+1)//2      # - (A+B) * t * S1(m)
 + A*B*m)                  # + A*B*m
// 2

which is exactly

\[ \frac{t^2 S_2(m) - (A + B)\,t\,S_1(m) + A B m}{2}, \]

where \[ S_1(m) = 1 + 2 + \dots + m = \frac{m(m+1)}{2}, \qquad S_2(m) = 1^2 + 2^2 + \dots + m^2 = \frac{m(m+1)(2m+1)}{6}. \]

This is the closed form for that double sum of floor functions derived from Legendre’s formula; the result is \( E_p \), the exponent of prime \( p \) in \( \mathrm{sdf}(N) \).

After the loop, exps[i] is the exponent \( E_{p_i} \) for primes[i].

So at this point we have, conceptually:

\[ \mathrm{sdf}(N) = \prod_i \text{primes}[i]^{\text{exps}[i]}. \]

Building an \(N\)-th root of unity

omega = pow(ROOT, (MOD - 1) // N, MOD)

Since ROOT is a primitive root modulo MOD, it has order \( \text{MOD} - 1 \).

\( \text{MOD} - 1 \) is divisible by \( N = 1000 \), so

\[ \omega = \text{ROOT}^{(\text{MOD} - 1)/N} \pmod{\text{MOD}} \]

is an element of order exactly \( N \).

That means \( \omega^N \equiv 1 \) and \( \omega^k \not\equiv 1 \) for \( 0 < k < N \), i.e. \( \omega \) is an \(N\)-th primitive root of unity in this finite field.

This is what we need for the root-of-unity filter.

Applying the root-of-unity filter

ans, w = 0, 1
for _ in range(N):
    prod = 1
    for p, e in zip(primes, exps):
        r = p * w % MOD
        term = (e + 1) % MOD if r == 1 else (pow(r, e + 1, MOD) - 1) * pow(r - 1, MOD - 2, MOD) % MOD
        prod = prod * term % MOD
    ans = (ans + prod) % MOD
    w = w * omega % MOD

This is computing

\[ \sum_{j=0}^{N-1} \prod_{p} \sum_{e=0}^{E_p} (p \omega^{j})^{e} \]

in modular arithmetic.

Details:

• The outer loop runs \( N \) times, with \( w = \omega^{j} \) at iteration \( j \): initially \( w = 1 = \omega^{0} \), and each step sets \( w \leftarrow w \cdot \omega \bmod \text{MOD} \), so \( w = \omega^{j+1} \).
• For each prime \( p \) and its exponent \( e = E_p \):
  • \( r = p \cdot w \bmod \text{MOD} = p \omega^j \).
  • We want the geometric sum \[ \sum_{k=0}^{E_p} r^k = \begin{cases} E_p + 1, & \text{if } r = 1,\\[4pt] \dfrac{r^{E_p+1} - 1}{r - 1}, & \text{if } r \ne 1. \end{cases} \]

That’s why the code does:

if r == 1:
    term = (e + 1) % MOD
else:
    term = (pow(r, e + 1, MOD) - 1) * pow(r - 1, MOD - 2, MOD) % MOD

Here pow(r, e + 1, MOD) is \( r^{E_p+1} \), and pow(r - 1, MOD - 2, MOD) is the modular inverse of \( r - 1 \) (by Fermat’s little theorem).

Then prod multiplies all these terms:

prod = prod * term % MOD

so after the inner loop:

\[ \text{prod} = \prod_{p} \sum_{k=0}^{E_p} (p \omega^{j})^{k}. \]

Then ans accumulates the sum over all \( j \):

ans = (ans + prod) % MOD

So when the loop finishes:

\[ \text{ans} \equiv \sum_{j=0}^{N-1} \prod_{p} \sum_{k=0}^{E_p} (p \omega^{j})^{k} \pmod{\text{MOD}}. \]

By the root-of-unity filter identity, this satisfies

\[ \text{ans} \equiv N \cdot D_0 \pmod{\text{MOD}}, \]

where \( D_0 \) is the desired “sum of divisors with total exponent sum \( \equiv 0 \pmod{N} \).”

Divide by N to finish

print(ans * pow(N, MOD - 2, MOD) % MOD)

Here pow(N, MOD - 2, MOD) is the modular inverse of \( N \) modulo \( \text{MOD} \).

So this computes

\[ D_0 \equiv \text{ans} \cdot N^{-1} \pmod{\text{MOD}}, \]

which is exactly the final answer required by the problem.

Last Word

• For N=5000: 975129331
• For N=10000: 749615570
• For N=25000: 872445465