Detailed mathematical approach

Problem Summary

Let \(\omega(n)\) denote the number of distinct prime divisors of \(n\). The problem uses the arithmetic function

$S(n)=\sum_{d\mid n} 2^{\omega(d)}.$

We must compute

$F(N)=\sum_{i=2}^{N} S(i!) \pmod{M},\qquad M=1{,}000{,}000{,}087,$

for the very large value \(N=10^7\). Directly factoring every factorial or enumerating divisors is far too slow, so the key is to update \(S(i!)\) from \(S((i-1)!)\) using only the prime factors of \(i\).

Mathematical Approach

The implementations rely on a multiplicative description of \(S(n)\), then apply it incrementally to the factorial sequence.

Step 1: Turn the divisor sum into a product

Write

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

Every divisor of \(n\) has the form

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

The value \(\omega(d)\) counts how many exponents \(a_p\) are positive. Therefore the divisor sum factors prime by prime. For one prime power \(p^e\), the local contribution is

$1+\underbrace{2+\cdots+2}_{e\text{ times}}=2e+1,$

because exponent \(0\) contributes \(2^{\omega(1)}=1\), while each exponent \(1,\dots,e\) contributes \(2^{\omega(p^a)}=2\).

Hence

$S(n)=\prod_{p^{e_p}\parallel n} (2e_p+1).$

As a quick check, \(6=2^1\cdot 3^1\), so

$S(6)=(2\cdot 1+1)(2\cdot 1+1)=3\cdot 3=9,$

which matches the divisor sum \(1+2+2+4=9\) over \(d=1,2,3,6\).

Step 2: Apply the product formula to factorials

For factorials, let

$E_p(i)=v_p(i!),$

so that

$i!=\prod_{p} p^{E_p(i)}.$

Then the previous identity becomes

$S(i!)=\prod_{p} \bigl(2E_p(i)+1\bigr).$

Now compare consecutive factorials:

$i!=i\cdot (i-1)!.$

If \(i\) contains prime \(p\) with multiplicity \(c=v_p(i)\), then

$E_p(i)=E_p(i-1)+c.$

All primes not dividing \(i\) keep the same exponent. This means only a few local factors change from one step to the next.

Step 3: Derive the incremental update formula

Suppose

$i=\prod_{p^{c_p}\parallel i} p^{c_p}.$

Each affected prime replaces one factor in the product for \(S((i-1)!)\):

$2E_p(i-1)+1 \longrightarrow 2\bigl(E_p(i-1)+c_p\bigr)+1.$

Therefore

$S(i!)=S((i-1)!)\prod_{p^{c_p}\parallel i}\frac{2\bigl(E_p(i-1)+c_p\bigr)+1}{2E_p(i-1)+1}.$

This is the central idea of the algorithm. It never recomputes the full product from scratch. It only divides out the old odd factor for each prime in \(i\), multiplies by the new odd factor, and keeps a running value of \(S(i!)\).

Step 4: Bound the needed odd factors and modular inverses

Every update uses odd numbers of the form

$2E_p(i)+1.$

The largest exponent in a factorial always belongs to prime \(2\), so for every prime \(p\) and every \(i\le N\),

$E_p(i)\le E_2(N)=v_2(N!).$

Hence all odd factors that can ever appear are at most

$B=2v_2(N!)+1.$

The implementation precomputes modular inverses up to this bound once, then reuses them during every ratio update. That turns each division in the formula above into a constant-time modular multiplication.

Step 5: Worked example

Take \(N=5\). The factorial values evolve as follows.

For \(2!=2\), we have \(2^1\), so

$S(2!)=2\cdot 1+1=3.$

For \(3!=6=2^1\cdot 3^1\),

$S(3!)=(2\cdot 1+1)(2\cdot 1+1)=3\cdot 3=9.$

Now move from \(3!\) to \(4!\). Since \(4=2^2\), the exponent of \(2\) rises from \(1\) to \(3\). The local factor for prime \(2\) changes from \(3\) to \(7\), so

$S(4!)=S(3!)\cdot \frac{7}{3}=9\cdot \frac{7}{3}=21.$

Next, \(5=5^1\), so prime \(5\) appears with exponent \(1\) for the first time. Its local factor changes from \(1\) to \(3\), giving

$S(5!)=S(4!)\cdot 3=21\cdot 3=63.$

Therefore

$F(5)=S(2!)+S(3!)+S(4!)+S(5!)=3+9+21+63=96.$

This example shows exactly why an incremental product update is much cheaper than refactoring each factorial independently.

How the Code Works

The C++, Python, and Java implementations first build a linear smallest-prime-factor table up to \(N\). That table lets the algorithm factor every integer \(i\) by repeated lookups and divisions, without trial division.

They also maintain an array of current factorial exponents \(E_p(i!)\) for all primes \(p\), together with one running modular value equal to the current \(S(i!)\). Before the main loop starts, the implementations compute the inverse table up to \(2v_2(N!)+1\), which covers every odd factor \(2E_p+1\) that may appear later.

For each \(i\) from \(2\) to \(N\), the factorization of \(i\) yields the multiplicity \(c_p\) of each prime dividing \(i\). The implementation updates the stored exponent of that prime, replaces the old local factor \(2e+1\) by the new factor \(2(e+c_p)+1\) inside the running product, and then adds the updated \(S(i!)\) to the cumulative total. No divisor list and no explicit factorial value is ever constructed.

Complexity Analysis

The linear smallest-prime-factor sieve costs \(O(N)\) time and \(O(N)\) memory. The inverse table also has size \(O(N)\), because \(v_2(N!)<N\).

Factoring every \(i\le N\) through the smallest-prime-factor table takes total time proportional to the total number of prime factors with multiplicity encountered along the way. On average this is \(O(N\log\log N)\), with a coarser upper bound of \(O(N\log N)\). The overall memory usage remains \(O(N)\).

Footnotes and References

1. Problem page: https://projecteuler.net/problem=675
2. Factorial: Wikipedia - Factorial
3. Prime factorization: Wikipedia - Prime factorization
4. Legendre's formula: Wikipedia - Legendre's formula
5. Multiplicative function: Wikipedia - Multiplicative function
6. Linear sieve: cp-algorithms - Linear Sieve