Problem 516: $5$-smooth Totients

Project Euler Problem 516 Solution

EulerSolve provides an optimized solution for Project Euler Problem 516, $5$-smooth Totients, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For a given limit $L$, define $S(L)$ as the sum of all integers $n\le L$ such that Euler's totient $\varphi(n)$ is 5-smooth, meaning that every prime factor of $\varphi(n)$ belongs to $\{2,3,5\}$. The target value is $S(10^{12}) \bmod 2^{32}$, so a direct scan of all $n\le 10^{12}$ is completely infeasible. Mathematical Approach The key observation is that the condition on $\varphi(n)$ forces a very rigid prime factorization for $n$. Once that structure is identified, the sum can be reorganized into a much smaller search over 5-smooth numbers and a recursively generated family of admissible prime products. Step 1: Separate the $2,3,5$ Part from the Larger Primes Write $$n=2^a3^b5^c\prod_{i=1}^r p_i^{e_i},\qquad p_i>5.$$ Euler's product formula gives $$\varphi(n)=\varphi(2^a)\varphi(3^b)\varphi(5^c)\prod_{i=1}^r p_i^{e_i-1}(p_i-1).$$ The factors coming from $2,3,5$ are always 5-smooth: $$\varphi(2^a)\in\{1,2,4,8,\dots\},\qquad \varphi(3^b)\in\{1,2,6,18,\dots\},\qquad \varphi(5^c)\in\{1,4,20,100,\dots\}.$$ So every possible obstruction comes from primes $p_i>5$. Step 2: Characterize Which Large Primes Are Allowed If some prime $p_i>5$ appears with exponent $e_i\ge 2$, then the factor $p_i^{e_i-1}$ divides $\varphi(n)$. That would force a prime factor larger than $5$ into $\varphi(n)$, which is impossible....

Detailed mathematical approach

Problem Summary

For a given limit $L$, define $S(L)$ as the sum of all integers $n\le L$ such that Euler's totient $\varphi(n)$ is 5-smooth, meaning that every prime factor of $\varphi(n)$ belongs to $\{2,3,5\}$. The target value is $S(10^{12}) \bmod 2^{32}$, so a direct scan of all $n\le 10^{12}$ is completely infeasible.

Mathematical Approach

The key observation is that the condition on $\varphi(n)$ forces a very rigid prime factorization for $n$. Once that structure is identified, the sum can be reorganized into a much smaller search over 5-smooth numbers and a recursively generated family of admissible prime products.

Step 1: Separate the $2,3,5$ Part from the Larger Primes

Write

$$n=2^a3^b5^c\prod_{i=1}^r p_i^{e_i},\qquad p_i>5.$$

Euler's product formula gives

$$\varphi(n)=\varphi(2^a)\varphi(3^b)\varphi(5^c)\prod_{i=1}^r p_i^{e_i-1}(p_i-1).$$

The factors coming from $2,3,5$ are always 5-smooth:

$$\varphi(2^a)\in\{1,2,4,8,\dots\},\qquad \varphi(3^b)\in\{1,2,6,18,\dots\},\qquad \varphi(5^c)\in\{1,4,20,100,\dots\}.$$

So every possible obstruction comes from primes $p_i>5$.

Step 2: Characterize Which Large Primes Are Allowed

If some prime $p_i>5$ appears with exponent $e_i\ge 2$, then the factor $p_i^{e_i-1}$ divides $\varphi(n)$. That would force a prime factor larger than $5$ into $\varphi(n)$, which is impossible. Therefore every prime $p>5$ may appear in $n$ at most once.

Moreover, when $p>5$ does appear, the factor $p-1$ divides $\varphi(n)$. Hence $p-1$ itself must be 5-smooth. This proves that every valid number has the shape

$$n=hq,$$

where $h=2^a3^b5^c$ is a 5-smooth number and $q$ is a squarefree product of distinct primes $p>5$ satisfying

$$p-1=2^\alpha 3^\beta 5^\gamma,\qquad \alpha,\beta,\gamma\ge 0.$$

The converse is also true. If $h$ is 5-smooth and $q=\prod p$ is a squarefree product of distinct primes with $p-1$ 5-smooth, then $h$ and $q$ are coprime, so

$$\varphi(hq)=\varphi(h)\prod_{p\mid q}(p-1),$$

which is again 5-smooth. Thus the decomposition is exact.

Step 3: Turn the Problem into a Sum over Admissible Prime Products

Let

$$\mathcal{H}(x)=\left\{2^a3^b5^c\le x:\ a,b,c\ge 0\right\}$$

be the set of 5-smooth numbers up to $x$, and define the prefix-sum function

$$A(x)=\sum_{h\in\mathcal{H}(x)} h.$$

Also define the admissible primes

$$\mathcal{P}(L)=\left\{p\le L:\ p>5,\ p\text{ prime},\ p-1\in\mathcal{H}(L)\right\}.$$

Every valid $n\le L$ is uniquely obtained by choosing a squarefree product $q$ of distinct primes from $\mathcal{P}(L)$, then choosing a 5-smooth multiplier $h\le L/q$. Therefore

$$S(L)=\sum_{q\in\mathcal{Q}(L)} \sum_{\substack{h\in\mathcal{H}(L)\\ h\le L/q}} qh,$$

where $\mathcal{Q}(L)$ is the set of squarefree products $q\le L$ built from distinct primes in $\mathcal{P}(L)$, including the empty product $q=1$.

Pulling $q$ outside the inner sum gives the main formula

$$\boxed{S(L)=\sum_{q\in\mathcal{Q}(L)} q\,A\left(\left\lfloor\frac{L}{q}\right\rfloor\right)\pmod{2^{32}}.}$$

Step 4: Why Enumerating $h+1$ Is Enough

The previous step shows that admissible primes are exactly the primes of the form

$$p=h+1,\qquad h\in\mathcal{H}(L).$$

So instead of searching through all primes up to $L$, we only generate 5-smooth numbers and test $h+1$ for primality. This is the decisive reduction: the 5-smooth list is tiny compared with the interval $[1,L]$.

After sorting the 5-smooth numbers, the values $A(x)$ can be answered by binary search and prefix sums. The remaining task is then to enumerate all feasible subset products $q\in\mathcal{Q}(L)$, which is done naturally by depth-first search with the pruning rule $q p\le L$.

Worked Example: $S(100)=3728$

For $L=100$, the 5-smooth numbers are

$$\{1,2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,27,30,32,36,40,45,48,50,54,60,64,72,75,80,81,90,96,100\}.$$

The admissible primes $p=h+1$ are

$$\{7,11,13,17,19,31,37,41,61,73,97\}.$$

Among their squarefree products, the only ones not exceeding $100$ are

$$1,\ 7,\ 11,\ 13,\ 17,\ 19,\ 31,\ 37,\ 41,\ 61,\ 73,\ 97,\ 77,\ 91.$$

The needed prefix sums are

$$A(100)=1258,\quad A(14)=60,\quad A(9)=38,\quad A(7)=21,\quad A(5)=15,\quad A(3)=6,\quad A(2)=3,\quad A(1)=1.$$

So the contributions are

$$\begin{aligned} q=1&:&&1\cdot A(100)=1258,\\ q=7&:&&7\cdot A(14)=420,\\ q=11&:&&11\cdot A(9)=418,\\ q=13&:&&13\cdot A(7)=273,\\ q=17&:&&17\cdot A(5)=255,\\ q=19&:&&19\cdot A(5)=285,\\ q=31&:&&31\cdot A(3)=186,\\ q=37&:&&37\cdot A(2)=111,\\ q=41&:&&41\cdot A(2)=123,\\ q=61&:&&61\cdot A(1)=61,\\ q=73&:&&73\cdot A(1)=73,\\ q=97&:&&97\cdot A(1)=97,\\ q=77&:&&77\cdot A(1)=77,\\ q=91&:&&91\cdot A(1)=91. \end{aligned}$$

Adding them gives

$$1258+420+418+273+255+285+186+111+123+61+73+97+77+91=3728,$$

which matches the checkpoint used by the implementation.

How the Code Works

The C++, Python, and Java implementations all use the same algorithm.

First, they generate every number of the form $2^a3^b5^c\le L$ by nested multiplicative loops. The resulting list is sorted, duplicates are removed, and a prefix-sum array modulo $2^{32}$ is built so that $A(x)$ can be queried quickly.

Next, for each 5-smooth number $h$, the implementation tests whether $h+1$ is prime. It uses fast modular exponentiation and a Miller-Rabin primality test with fixed small bases, which is sufficient for the numeric range of this problem. Every prime $h+1>5$ is stored as an admissible prime.

Finally, a depth-first search enumerates subset products $q$ of admissible primes in increasing order. At each recursive state it adds

$$q\,A\left(\left\lfloor\frac{L}{q}\right\rfloor\right)\pmod{2^{32}}$$

to the answer, then tries to append later admissible primes as long as the product stays at most $L$. Because the search only moves forward through the sorted prime list, each squarefree subset is visited exactly once.

Complexity Analysis

Let $H=|\mathcal{H}(L)|$, let $P=|\mathcal{P}(L)|$, and let $T=|\mathcal{Q}(L)|$, the number of subset products actually visited by the depth-first search. Generating the 5-smooth list takes $O(H)$ multiplicative steps, while sorting and deduplicating costs $O(H\log H)$. Building prefix sums is linear in $H$.

Primality is tested only on the $H$ candidates $h+1$, not on every integer up to $L$. With a fixed-base Miller-Rabin test, this contributes roughly $O(H\log L)$ modular-arithmetic work. The depth-first search visits each feasible subset product once and performs one binary search on the 5-smooth list per node, so the summation phase costs $O(T\log H)$.

Thus the algorithm is dominated by the sizes of the 5-smooth list and the feasible subset-product tree, both of which are tiny compared with $L$ itself. Memory usage is $O(H+P)$, plus recursion depth proportional to the number of primes currently selected.

Footnotes and References

Problem page: https://projecteuler.net/problem=516
Euler's totient function: Wikipedia — Euler's totient function
Regular numbers / 5-smooth numbers: Wikipedia — Regular number
Miller-Rabin primality test: Wikipedia — Miller-Rabin primality test
Squarefree integers: Wikipedia — Square-free integer

Problem 516 source code

C++

#include <algorithm>
#include <cstdint>
#include <functional>
#include <iostream>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = __uint128_t;

constexpr u64 kMask32 = 0xFFFF'FFFFULL;

u64 mod_mul(const u64 a, const u64 b, const u64 mod) {
    return static_cast<u64>((static_cast<u128>(a) * b) % mod);
}

u64 mod_pow(u64 base, u64 exp, const u64 mod) {
    u64 result = 1ULL % mod;
    u64 cur = base % mod;
    u64 e = exp;
    while (e > 0ULL) {
        if (e & 1ULL) {
            result = mod_mul(result, cur, mod);
        }
        cur = mod_mul(cur, cur, mod);
        e >>= 1ULL;
    }
    return result;
}

bool is_prime(const u64 n) {
    if (n < 2ULL) {
        return false;
    }
    for (const u64 p : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL}) {
        if (n == p) {
            return true;
        }
        if (n % p == 0ULL) {
            return false;
        }
    }

    u64 d = n - 1ULL;
    int s = 0;
    while ((d & 1ULL) == 0ULL) {
        d >>= 1ULL;
        ++s;
    }

    for (const u64 a : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL}) {
        if (a >= n) {
            continue;
        }
        u64 x = mod_pow(a, d, n);
        if (x == 1ULL || x == n - 1ULL) {
            continue;
        }
        bool witness = true;
        for (int r = 1; r < s; ++r) {
            x = mod_mul(x, x, n);
            if (x == n - 1ULL) {
                witness = false;
                break;
            }
        }
        if (witness) {
            return false;
        }
    }
    return true;
}

std::vector<u64> generate_hamming(const u64 limit) {
    std::vector<u64> out;
    for (u64 a = 1ULL; a <= limit; a *= 2ULL) {
        for (u64 b = a; b <= limit; b *= 3ULL) {
            for (u64 c = b; c <= limit; c *= 5ULL) {
                out.push_back(c);
                if (c > limit / 5ULL) {
                    break;
                }
            }
            if (b > limit / 3ULL) {
                break;
            }
        }
        if (a > limit / 2ULL) {
            break;
        }
    }
    std::sort(out.begin(), out.end());
    out.erase(std::unique(out.begin(), out.end()), out.end());
    return out;
}

u64 solve(const u64 limit) {
    const std::vector<u64> hamming = generate_hamming(limit);

    std::vector<u64> prefix_mod(static_cast<std::size_t>(hamming.size() + 1), 0ULL);
    for (std::size_t i = 0; i < hamming.size(); ++i) {
        prefix_mod[i + 1] = (prefix_mod[i] + (hamming[i] & kMask32)) & kMask32;
    }

    std::vector<u64> admissible_primes;
    admissible_primes.reserve(hamming.size());
    for (const u64 h : hamming) {
        const u64 p = h + 1ULL;
        if (p > 5ULL && p <= limit && is_prime(p)) {
            admissible_primes.push_back(p);
        }
    }
    std::sort(admissible_primes.begin(), admissible_primes.end());

    u64 total_mod = 0ULL;

    const auto sum_hamming_mod = [&](const u64 x) -> u64 {
        const auto it = std::upper_bound(hamming.begin(), hamming.end(), x);
        const std::size_t count = static_cast<std::size_t>(it - hamming.begin());
        return prefix_mod[count];
    };

    std::function<void(std::size_t, u64)> dfs = [&](const std::size_t idx, const u64 prod) {
        const u64 h_sum = sum_hamming_mod(limit / prod);
        total_mod = (total_mod + ((prod & kMask32) * h_sum & kMask32)) & kMask32;

        for (std::size_t i = idx; i < admissible_primes.size(); ++i) {
            const u64 p = admissible_primes[i];
            if (prod > limit / p) {
                break;
            }
            dfs(i + 1, prod * p);
        }
    };

    dfs(0, 1ULL);
    return total_mod;
}

u64 brute(const int limit) {
    std::vector<int> phi(static_cast<std::size_t>(limit + 1), 0);
    for (int i = 0; i <= limit; ++i) {
        phi[static_cast<std::size_t>(i)] = i;
    }
    for (int p = 2; p <= limit; ++p) {
        if (phi[static_cast<std::size_t>(p)] != p) {
            continue;
        }
        for (int m = p; m <= limit; m += p) {
            phi[static_cast<std::size_t>(m)] -= phi[static_cast<std::size_t>(m)] / p;
        }
    }

    auto is_smooth = [](int x) -> bool {
        while (x % 2 == 0) {
            x /= 2;
        }
        while (x % 3 == 0) {
            x /= 3;
        }
        while (x % 5 == 0) {
            x /= 5;
        }
        return x == 1;
    };

    u64 sum_mod = 0ULL;
    for (int n = 1; n <= limit; ++n) {
        if (is_smooth(phi[static_cast<std::size_t>(n)])) {
            sum_mod = (sum_mod + static_cast<u64>(n)) & kMask32;
        }
    }
    return sum_mod;
}

bool run_checkpoints() {
    if (solve(100ULL) != 3'728ULL) {
        std::cerr << "Checkpoint failed: S(100)\n";
        return false;
    }
    if (solve(5'000ULL) != brute(5'000)) {
        std::cerr << "Checkpoint failed: formula/bruteforce mismatch at 5000\n";
        return false;
    }
    return true;
}

}  // namespace

int main() {
    if (!run_checkpoints()) {
        return 1;
    }

    constexpr u64 limit = 1'000'000'000'000ULL;
    std::cout << solve(limit) << '\n';
    return 0;
}

Python

import math

def solve():
    LIMIT = 10**12
    MASK32 = 0xFFFFFFFF

    def mod_mul(a, b, mod):
        return a * b % mod

    def mod_pow(base, exp, mod):
        result = 1 % mod
        cur = base % mod
        while exp > 0:
            if exp & 1: result = mod_mul(result, cur, mod)
            cur = mod_mul(cur, cur, mod)
            exp >>= 1
        return result

    def is_prime(n):
        if n < 2: return False
        for p in [2, 3, 5, 7, 11, 13]:
            if n == p: return True
            if n % p == 0: return False
        d = n - 1
        s = 0
        while d % 2 == 0: d >>= 1; s += 1
        for a in [2, 3, 5, 7, 11, 13]:
            if a >= n: continue
            x = mod_pow(a, d, n)
            if x == 1 or x == n - 1: continue
            witness = True
            for _ in range(1, s):
                x = mod_mul(x, x, n)
                if x == n - 1: witness = False; break
            if witness: return False
        return True

    def gen_hamming(limit):
        out = []
        a = 1
        while a <= limit:
            b = a
            while b <= limit:
                c = b
                while c <= limit:
                    out.append(c)
                    if c > limit // 5: break
                    c *= 5
                if b > limit // 3: break
                b *= 3
            if a > limit // 2: break
            a *= 2
        out = sorted(set(out))
        return out

    hamming = gen_hamming(LIMIT)
    prefix_mod = [0] * (len(hamming) + 1)
    for i in range(len(hamming)):
        prefix_mod[i+1] = (prefix_mod[i] + (hamming[i] & MASK32)) & MASK32

    admissible = sorted(p for h in hamming if (p := h + 1) > 5 and p <= LIMIT and is_prime(p))

    import bisect
    def sum_hamming_mod(x):
        cnt = bisect.bisect_right(hamming, x)
        return prefix_mod[cnt]

    total_mod = 0
    def dfs(idx, prod):
        nonlocal total_mod
        h_sum = sum_hamming_mod(LIMIT // prod)
        total_mod = (total_mod + ((prod & MASK32) * h_sum & MASK32)) & MASK32
        for i in range(idx, len(admissible)):
            p = admissible[i]
            if prod > LIMIT // p: break
            dfs(i + 1, prod * p)

    dfs(0, 1)
    return str(total_mod)

if __name__ == '__main__':
    print(solve())

Java

import java.util.ArrayList;
import java.util.Collections;
import java.util.List;

public class Euler516 {

    static long modPow(long base, long exp, long mod) {
        long result = 1 % mod;
        long cur = base % mod;
        long e = exp;
        while (e > 0) {
            if ((e & 1) == 1) {
                result = modMul(result, cur, mod);
            }
            cur = modMul(cur, cur, mod);
            e >>= 1;
        }
        return result;
    }

    static long modMul(long a, long b, long mod) {
        long q = (long) ((double) a * b / mod);
        long r = a * b - q * mod;
        while (r < 0)
            r += mod;
        while (r >= mod)
            r -= mod;
        return r;
    }

    static boolean isPrime(long n) {
        if (n < 2)
            return false;
        long[] p_list = { 2, 3, 5, 7, 11, 13 };
        for (long p : p_list) {
            if (n == p)
                return true;
            if (n % p == 0)
                return false;
        }

        long d = n - 1;
        int s = 0;
        while ((d & 1) == 0) {
            d >>= 1;
            s++;
        }

        for (long a : p_list) {
            if (a >= n)
                break;
            long x = modPow(a, d, n);
            if (x == 1 || x == n - 1)
                continue;
            boolean witness = true;
            for (int r = 1; r < s; r++) {
                x = modMul(x, x, n);
                if (x == n - 1) {
                    witness = false;
                    break;
                }
            }
            if (witness)
                return false;
        }
        return true;
    }

    static List<Long> generateHamming(long limit) {
        List<Long> out = new ArrayList<>();
        for (long a = 1; a <= limit; a *= 2) {
            for (long b = a; b <= limit; b *= 3) {
                for (long c = b; c <= limit; c *= 5) {
                    out.add(c);
                    if (c > limit / 5)
                        break;
                }
                if (b > limit / 3)
                    break;
            }
            if (a > limit / 2)
                break;
        }
        Collections.sort(out);
        List<Long> uniqueOut = new ArrayList<>();
        if (!out.isEmpty()) {
            uniqueOut.add(out.get(0));
            for (int i = 1; i < out.size(); i++) {
                if (!out.get(i).equals(out.get(i - 1))) {
                    uniqueOut.add(out.get(i));
                }
            }
        }
        return uniqueOut;
    }

    static final long kMask32 = 0xFFFFFFFFL;
    static long totalMod = 0;
    static long limit = 1000000000000L;
    static List<Long> hamming;
    static long[] prefixMod;
    static List<Long> admissiblePrimes;

    static long sumHammingMod(long x) {
        int low = 0, high = hamming.size();
        while (low < high) {
            int mid = low + (high - low) / 2;
            if (hamming.get(mid) <= x) {
                low = mid + 1;
            } else {
                high = mid;
            }
        }
        return prefixMod[low];
    }

    static void dfs(int idx, long prod) {
        long hSum = sumHammingMod(limit / prod);
        totalMod = (totalMod + ((prod & kMask32) * hSum & kMask32)) & kMask32;

        for (int i = idx; i < admissiblePrimes.size(); i++) {
            long p = admissiblePrimes.get(i);
            if (prod > limit / p)
                break;
            dfs(i + 1, prod * p);
        }
    }

    public static void main(String[] args) {
        hamming = generateHamming(limit);
        prefixMod = new long[hamming.size() + 1];
        for (int i = 0; i < hamming.size(); i++) {
            prefixMod[i + 1] = (prefixMod[i] + (hamming.get(i) & kMask32)) & kMask32;
        }

        admissiblePrimes = new ArrayList<>();
        for (long h : hamming) {
            long p = h + 1;
            if (p > 5 && p <= limit && isPrime(p)) {
                admissiblePrimes.add(p);
            }
        }
        Collections.sort(admissiblePrimes);

        dfs(0, 1);

        System.out.println(totalMod);
    }
}