Solution

Answer: 476001479068717

For a positive integer $n$, define $s(n)$ as the smallest integer $m$ such that

$$n \mid m!.$$

We must compute

$$S(10^8)=\sum_{n=2}^{10^8} s(n).$$

The key is to understand $s(n)$ through prime factorizations.

1. Mathematical analysis

Suppose

$$n=\prod p_i^{a_i}.$$

Then $n\mid m!$ iff for every prime power $p_i^{a_i}$,

$$v_{p_i}(m!)\ge a_i,$$

where $v_p(x)$ is the exponent of $p$ in $x$.

Thus:

$$s(n)=\max_i s(p_i^{a_i}),$$

so the entire problem reduces to prime powers.

Prime powers

For a fixed prime $p$, define:

$$f(p,a)=s(p^a).$$

We need the smallest $m$ such that

$$v_p(m!)\ge a.$$

By Legendre’s formula:

$v_p(m!)=\sum_{k\ge1}\left\lfloor\frac{m}{p^k}\right\rfloor$

So $f(p,a)$ is the minimal $m$ satisfying this inequality.

Examples:

$f(5,1)=5$
$f(5,2)=10$ because $10!$ contains exactly two factors of $5$

Computing all $s(n)$

For every prime power $q=p^a\le N$:

Compute $m=f(p,a)$
Every multiple of $q$ must satisfy $s(k)\ge m$

Therefore we can sieve:

$$s(k)=\max(s(k),,f(p,a)) \quad\text{for all }k\equiv0\pmod{p^a}.$$

This is very similar to a prime sieve.

The total complexity is approximately

$$N\sum_{p^a\le N}\frac1{p^a},$$

which is fast enough for $N=10^8$ in optimized code.

2. Python implementation

from math import isqrt

N = 10**8

# ------------------------------------------------------------
# Legendre valuation: exponent of prime p inside m!
# ------------------------------------------------------------
def vp_fact(m, p):
    s = 0
    while m:
        m //= p
        s += m
    return s

# ------------------------------------------------------------
# Smallest m such that v_p(m!) >= a
# ------------------------------------------------------------
def min_m(p, a):
    lo, hi = 1, p * a

    while lo < hi:
        mid = (lo + hi) // 2

        if vp_fact(mid, p) >= a:
            hi = mid
        else:
            lo = mid + 1

    return lo

# ------------------------------------------------------------
# Sieve primes up to N
# ------------------------------------------------------------
is_comp = bytearray(N + 1)
primes = []

for i in range(2, N + 1):
    if not is_comp[i]:
        primes.append(i)

        if i <= isqrt(N):
            step = i
            start = i * i
            for j in range(start, N + 1, step):
                is_comp[j] = 1

# ------------------------------------------------------------
# Compute s(n) for all n <= N
# ------------------------------------------------------------
s = [0] * (N + 1)

for p in primes:

    power = p
    a = 1

    while power <= N:

        m = min_m(p, a)

        # every multiple of p^a needs at least m
        for k in range(power, N + 1, power):
            if s[k] < m:
                s[k] = m

        if power > N // p:
            break

        power *= p
        a += 1

# ------------------------------------------------------------
# Sum all values
# ------------------------------------------------------------
answer = sum(s[2:])
print(answer)

3. Code walkthrough

`vp_fact`

def vp_fact(m, p):

Implements Legendre’s formula:

$$v_p(m!)=\left\lfloor\frac{m}{p}\right\rfloor +\left\lfloor\frac{m}{p^2}\right\rfloor+\cdots$$

`min_m`

Binary searches the smallest $m$ with enough powers of $p$.

Example:

For $p=5,a=2$:
$v_5(9!)=1$
$v_5(10!)=2$

So min_m(5,2)=10.

Prime sieve

Generates all primes up to $10^8$.

Main sieve

For every prime power $p^a$:

for k in range(power, N + 1, power):

every multiple of $p^a$ must have

$$s(k)\ge s(p^a).$$

We keep the maximum over all prime powers dividing $k$.

That exactly equals:

$$s(k)=\max_{p^a\mid k} s(p^a).$$

4. Verification on the sample

For $N=100$, the program produces:

$$S(100)=2012,$$

matching the statement.

The final result is about $4.76\times10^{14}$, which is reasonable because average values of $s(n)$ are on the order of millions over $10^8$ terms.

Answer: 476001479068717

Project Euler Problem 549