Problem 516: 5-Smooth Totients

Mathematical Foundation

Theorem 1 (Characterization of 5-Smooth Totient Numbers). A positive integer $n$ satisfies ”$\varphi(n)$ is 5-smooth” if and only if

where $a, b, c \geq 0$, $k \geq 0$, and $p_1, p_2, \ldots, p_k$ are distinct primes (each $> 5$) such that $p_i - 1$ is a Hamming number for every $i$.

Proof. ($\Leftarrow$) By multiplicativity of $\varphi$ on pairwise coprime factors:

We have $\varphi(2^a) = 2^{a-1}$ (for $a \geq 1$; $\varphi(1) = 1$), $\varphi(3^b) = 2 \cdot 3^{b-1}$, $\varphi(5^c) = 4 \cdot 5^{c-1}$, and $\varphi(p_i) = p_i - 1$, each of which is 5-smooth by hypothesis. A product of 5-smooth numbers is 5-smooth, so $\varphi(n)$ is 5-smooth.

($\Rightarrow$) Suppose $\varphi(n)$ is 5-smooth. Let $q > 5$ be a prime dividing $n$.

• If $q^2 \mid n$, then $\varphi(q^2) = q(q-1)$ divides $\varphi(n)$. Since $q > 5$, the factor $q$ makes $\varphi(n)$ non-5-smooth, a contradiction. So $q \| n$ (exact power 1).
• Since $q \| n$, we have $\varphi(q) = q - 1 \mid \varphi(n)$. For $\varphi(n)$ to be 5-smooth, $q - 1$ must be 5-smooth.

Thus every prime $q > 5$ dividing $n$ appears to the first power and satisfies $q - 1$ being 5-smooth. The remaining part of $n$ is $2^a \cdot 3^b \cdot 5^c$, and $\varphi$ of any prime power of 2, 3, or 5 is 5-smooth. $\square$

Lemma 1 (Hamming Number Count). The number of Hamming numbers not exceeding $L$ is $O(\log^3 L)$.

Proof. A Hamming number $h = 2^a \cdot 3^b \cdot 5^c \leq L$ requires $a \leq \log_2 L$, $b \leq \log_3 L$, $c \leq \log_5 L$. The count is at most $(\log_2 L + 1)(\log_3 L + 1)(\log_5 L + 1) = O(\log^3 L)$. For $L = 10^{12}$, this is approximately 3,428. $\square$

Lemma 2 (Hamming Prime Count). A Hamming prime is a prime $p$ such that $p - 1$ is a Hamming number. The number of Hamming primes up to $L$ is at most the number of Hamming numbers up to $L$, i.e., $O(\log^3 L)$. For $L = 10^{12}$, there are approximately 545 such primes.

Proof. Each Hamming prime $p$ corresponds to a unique Hamming number $p - 1$, so the count is bounded by the number of Hamming numbers. $\square$

Editorial

We generate all Hamming numbers up to L. We then find Hamming primes. Finally, iterate over h in hamming. We perform a recursive search over the admissible choices, prune branches that violate the derived constraints, and keep only the candidates that satisfy the final condition.

Pseudocode

Generate all Hamming numbers up to L
Find Hamming primes
for h in hamming
Enumerate valid n via DFS over subsets of Hamming primes
Multiply product by each Hamming number h with product * h <= L
and add product * h to total
for h in hamming
Extend product by including more Hamming primes

Complexity Analysis

• Time: The Hamming number generation is $O(\log^3 L)$. Primality testing for each Hamming number costs $O(\sqrt{h})$ per test, totaling $O(\log^3 L \cdot \sqrt{L})$ in the worst case (mitigated by Miller-Rabin). The DFS over subsets of $\sim$545 primes is pruned aggressively by the $\leq L$ bound; in practice, the tree is small and the total runtime is under 1 second.
• Space: $O(\log^3 L)$ for the Hamming number list, plus $O(k)$ recursion depth where $k$ is the number of Hamming primes used (at most $\sim$40 for $L = 10^{12}$).

Answer

#include <bits/stdc++.h>
using namespace std;

typedef unsigned long long ull;
typedef long long ll;
const ull LIMIT = 1000000000000ULL; // 10^12
const ull MOD = (1ULL << 32);

// Miller-Rabin primality test for numbers up to 10^12
ull mulmod(ull a, ull b, ull m) {
    return (__uint128_t)a * b % m;
}

ull powmod(ull a, ull b, ull m) {
    ull res = 1;
    a %= m;
    while (b > 0) {
        if (b & 1) res = mulmod(res, a, m);
        a = mulmod(a, a, m);
        b >>= 1;
    }
    return res;
}

bool millerRabin(ull n, ull a) {
    if (n % a == 0) return n == a;
    ull d = n - 1;
    int r = 0;
    while (d % 2 == 0) { d /= 2; r++; }
    ull x = powmod(a, d, n);
    if (x == 1 || x == n - 1) return true;
    for (int i = 0; i < r - 1; i++) {
        x = mulmod(x, x, n);
        if (x == n - 1) return true;
    }
    return false;
}

bool isPrime(ull n) {
    if (n < 2) return false;
    if (n < 4) return true;
    if (n % 2 == 0 || n % 3 == 0) return false;
    // Deterministic for n < 3.3 * 10^24
    for (ull a : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
        if (!millerRabin(n, a)) return false;
    }
    return true;
}

int main() {
    // Step 1: Generate all Hamming numbers <= LIMIT
    vector<ull> hamming;
    for (ull a = 1; a <= LIMIT; a *= 2) {
        for (ull b = a; b <= LIMIT; b *= 3) {
            for (ull c = b; c <= LIMIT; c *= 5) {
                hamming.push_back(c);
            }
        }
    }
    sort(hamming.begin(), hamming.end());

    // Step 2: Find Hamming primes (h+1 is prime where h is a Hamming number)
    vector<ull> hprimes;
    for (ull h : hamming) {
        if (h >= 2 && isPrime(h + 1)) {
            hprimes.push_back(h + 1);
        }
    }
    // Add 2 (since phi(2) = 1 which is Hamming, but 2-1=1 is Hamming)
    // Actually 2 = 1+1, and 1 is a Hamming number, so 2 is already included
    // We also need to handle that 2,3,5 themselves are special
    // Remove 2, 3, 5 from hprimes since they are handled as part of Hamming numbers
    // Actually: n can have powers of 2,3,5 AND distinct Hamming primes
    // The primes 2,3,5 with higher powers are already covered by Hamming numbers
    // For primes p where p-1 is Hamming: p=2(1+1),3(2+1),5(4+1),7(6+1=2*3),11,13,17,...
    // We keep all but note that 2,3,5 are special (can appear to any power)

    // Remove 2, 3, 5 from hprimes - they're handled separately as Hamming base
    hprimes.erase(remove_if(hprimes.begin(), hprimes.end(),
        [](ull p) { return p == 2 || p == 3 || p == 5; }), hprimes.end());
    sort(hprimes.begin(), hprimes.end());

    // Step 3: For each subset of Hamming primes (product <= LIMIT),
    // multiply by each Hamming number <= LIMIT/product and sum up
    ull ans = 0;

    // DFS to enumerate subsets of hprimes
    // For each product of a subset, multiply by all Hamming numbers that fit
    function<void(int, ull)> dfs = [&](int idx, ull prod) {
        // Add contribution: prod * each Hamming number <= LIMIT/prod
        ull maxH = LIMIT / prod;
        for (ull h : hamming) {
            if (h > maxH) break;
            ans = (ans + (prod * h) % MOD) % MOD;
        }
        // Try adding the next prime
        for (int i = idx; i < (int)hprimes.size(); i++) {
            if (prod > LIMIT / hprimes[i]) break;
            dfs(i + 1, prod * hprimes[i]);
        }
    };

    dfs(0, 1);

    cout << ans << endl;
    return 0;
}

"""
Project Euler Problem 516: 5-Smooth Totients

Find S(10^12) mod 2^32, where S(L) is the sum of all n <= L
such that phi(n) is a 5-smooth (Hamming) number.
"""

from sympy import isprime

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

def solve():
    # Step 1: Generate all Hamming 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.sort()

    # Step 2: Find Hamming primes (h+1 is prime, h is Hamming, excluding 2,3,5)
    hprimes = []
    for h in hamming:
        if h >= 2 and isprime(h + 1):
            p = h + 1
            if p not in (2, 3, 5):
                hprimes.append(p)
    hprimes.sort()

    # Step 3: DFS over subsets of Hamming primes, multiply by Hamming numbers
    ans = 0

    def dfs(idx, prod):
        nonlocal ans
        max_h = LIMIT // prod
        for h in hamming:
            if h > max_h:
                break
            ans = (ans + prod * h) % MOD

        for i in range(idx, len(hprimes)):
            if prod > LIMIT // hprimes[i]:
                break
            dfs(i + 1, prod * hprimes[i])

    dfs(0, 1)
    print(ans)

if __name__ == "__main__":
    solve()