Project Euler Problem 694

Solution

Answer: 1339784153569958487

Let $C$ be the set of cube-full numbers. Since $s(n)$ counts cube-full divisors of $n$,

$$S(N)=\sum_{n\le N}s(n)$$

can be re-indexed by divisors:

$$S(N) = \sum_{d\in C,\ d\le N} \left\lfloor \frac{N}{d}\right\rfloor$$

because each cube-full $d$ contributes $1$ to every multiple of $d$.

So the task reduces to enumerating all cube-full numbers $d\le 10^{18}$ and summing $\lfloor N/d\rfloor$.

A cube-full number has every prime exponent at least $3$, so we can recursively generate them uniquely by multiplying prime powers $p^e$ with $e\ge 3$.

Python implementation used:

from math import isqrt

def sieve(n):
    bs = bytearray(b'\x01') * (n + 1)
    bs[:2] = b'\x00\x00'
    for i in range(2, isqrt(n) + 1):
        if bs[i]:
            bs[i * i:n + 1:i] = b'\x00' * (((n - i * i) // i) + 1)
    return [i for i, b in enumerate(bs) if b]

primes = sieve(10**6)  # cube root of 10^18

def S(N):
    total = N  # contribution from d = 1

    def dfs(start_index, current):
        nonlocal total
        remaining = N // current

        for j in range(start_index, len(primes)):
            p = primes[j]
            p3 = p * p * p
            if p3 > remaining:
                break

            pe = p3
            while pe <= remaining:
                v = current * pe
                total += N // v
                dfs(j + 1, v)
                pe *= p

    dfs(0, 1)
    return total

print(S(16))      # 19
print(S(100))     # 126
print(S(10000))   # 13344
print(S(10**18))

Verification:

$S(16)=19$ ✓
$S(100)=126$ ✓
$S(10000)=13344$ ✓

The computation for $10^{18}$ yields:

Answer: 1339784153569958487