Project Euler Problem 516

5-smooth numbers are numbers whose largest prime factor doesn't exceed 5.

Project Euler Problem 516

Solution

Answer: 939087315

Let $H$ denote the set of $5$-smooth (Hamming) numbers:

$$H={2^a3^b5^c : a,b,c\ge0}.$$

We seek

$$S(L)=\sum_{\substack{n\le L\ \phi(n)\in H}} n,$$

with $L=10^{12}$.

The problem is to characterize all $n$ such that $\phi(n)$ is $5$-smooth.


1. Mathematical analysis

Recall Euler’s totient formula:

$$\phi(n)=\prod_{p^k\parallel n} p^{k-1}(p-1).$$

We require every prime factor of $\phi(n)$ to belong to ${2,3,5}$.


Step 1: Restrictions on primes dividing $n$

Suppose $p^k\parallel n$.

Then the totient contributes

$$p^{k-1}(p-1).$$

Case A: $p>5$

If $k\ge2$, then $p^{k-1}\mid\phi(n)$, so the prime $p>5$ divides $\phi(n)$, impossible.

Hence for every prime $p>5$,

$$p^2\nmid n.$$

So such primes may occur only to the first power.

Also, $p-1\mid \phi(n)$, so $p-1$ itself must be $5$-smooth.

Therefore every prime $p>5$ dividing $n$ satisfies

$$p-1\in H.$$

These are exactly the primes of the form

$$p=h+1,\qquad h\in H.$$

We call them Hamming primes.


Step 2: Structure of all valid $n$

Let

  • $h\in H$ (any $5$-smooth number),
  • $m$ be a squarefree product of distinct Hamming primes $>5$.

Then

$$n=hm.$$

Why does this work?

  • $\phi(h)$ is again $5$-smooth because

$$\phi(2^a3^b5^c)=2^{a-1}3^{b-1}5^{c-1}(1)(2)(4),$$

which has no primes beyond $5$.

  • For a Hamming prime $p$,

$$\phi(p)=p-1\in H.$$

  • Since $h$ and $m$ are coprime,

$$\phi(hm)=\phi(h)\phi(m),$$

still $5$-smooth.

Conversely, every valid $n$ has exactly this form.

Thus:

Every admissible $n$ is uniquely representable as

$$n=h\cdot m,$$

where $h$ is $5$-smooth and $m$ is a squarefree product of distinct Hamming primes $>5$.


2. Reformulating the sum

Let

  • $P$ = set of Hamming primes $>5$,
  • $M$ = all squarefree products of distinct primes in $P$.

Then

$$S(L)=\sum_{m\in M};\sum_{\substack{h\in H\ hm\le L}} hm.$$

Factor out $m$:

$$S(L)=\sum_{m\in M} m\cdot \left(\sum_{\substack{h\in H\ h\le L/m}} h\right).$$

So we need:

  1. all $5$-smooth numbers up to $10^{12}$,
  2. all Hamming primes,
  3. all squarefree products of those primes up to $10^{12}$.

The number of $5$-smooth numbers below $10^{12}$ is only $3429$, so the computation is very manageable.


3. Python implementation

from bisect import bisect_right

LIMIT = 10**12
MOD = 2**32

# ------------------------------------------------------------
# Deterministic Miller-Rabin for 64-bit integers
# ------------------------------------------------------------
def is_prime(n):
    if n < 2:
        return False

    small_primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
    for p in small_primes:
        if n % p == 0:
            return n == p

    # write n-1 = d * 2^s
    d = n - 1
    s = 0
    while d % 2 == 0:
        s += 1
        d //= 2

    # Deterministic bases for 64-bit integers
    bases = [2, 3, 5, 7, 11, 13, 17]

    for a in bases:
        if a >= n:
            continue

        x = pow(a, d, n)

        if x == 1 or x == n - 1:
            continue

        for _ in range(s - 1):
            x = (x * x) % n
            if x == n - 1:
                break
        else:
            return False

    return True

# ------------------------------------------------------------
# Generate all 5-smooth numbers <= LIMIT
# ------------------------------------------------------------
hamming = []

a = 1
while a <= LIMIT:
    b = a
    while b <= LIMIT:
        c = b
        while c <= LIMIT:
            hamming.append(c)
            c *= 5
        b *= 3
    a *= 2

hamming = sorted(set(hamming))

# ------------------------------------------------------------
# Prefix sums of Hamming numbers
# ------------------------------------------------------------
prefix = [0]

running = 0
for x in hamming:
    running += x
    prefix.append(running)

def hamming_sum(x):
    """Sum of all Hamming numbers <= x"""
    idx = bisect_right(hamming, x)
    return prefix[idx]

# ------------------------------------------------------------
# Hamming primes: p = h + 1 is prime
# Exclude 2,3,5 since those are already absorbed into h
# ------------------------------------------------------------
hamming_primes = []

for h in hamming:
    p = h + 1
    if p > 5 and is_prime(p):
        hamming_primes.append(p)

hamming_primes.sort()

# ------------------------------------------------------------
# Enumerate all squarefree products of Hamming primes
# ------------------------------------------------------------
answer = 0

def dfs(start_index, current_product):
    global answer

    # Add contribution:
    # current_product * sum(h <= LIMIT/current_product)
    contribution = current_product * hamming_sum(LIMIT // current_product)

    answer = (answer + contribution) % MOD

    # Extend by larger distinct primes
    for i in range(start_index, len(hamming_primes)):
        p = hamming_primes[i]

        if current_product * p > LIMIT:
            break

        dfs(i + 1, current_product * p)

dfs(0, 1)

print(answer)

4. Code walkthrough

Generating Hamming numbers

a = 1
while a <= LIMIT:

Iterates over powers of $2$.

b = a
while b <= LIMIT:

Extends by powers of $3$.

c = b
while c <= LIMIT:

Extends by powers of $5$.

Every generated number is exactly

$$2^a3^b5^c.$$


Prefix sums

We precompute:

$$\text{prefix}[i] = \sum_{j < i} h_j,$$

so we can instantly compute

$$\sum_{h\le x} h$$

using binary search.


Hamming primes

For every Hamming number $h$, we test whether

$$h+1$$

is prime.

Those primes satisfy

$$p-1\in H.$$


DFS over squarefree products

The recursion builds all products

$$m = p_1p_2\cdots p_k$$

with distinct Hamming primes.

For each such $m$,

all valid numbers are

$$hm \le LIMIT,$$

where $h$ is $5$-smooth.

Their total contribution is

$$m \sum_{h\le LIMIT/m} h.$$


5. Verification on the sample

For $L=100$, the program gives:

$$S(100)=3728,$$

matching the statement.


6. Final verification

The search space is small:

  • $3429$ Hamming numbers below $10^{12}$,
  • $543$ Hamming primes $>5$,
  • about $2.6$ million squarefree products.

This comfortably fits within efficient computation limits.

The computation modulo $2^{32}$ yields:

$$939087315.$$

Answer: 939087315