Problem 204: Generalised Hamming Numbers

Project Euler Problem 204 Solution

EulerSolve provides an optimized solution for Project Euler Problem 204, Generalised Hamming Numbers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary A generalized Hamming number of type $T$ is a positive integer whose prime divisors are all at most $T$. For this problem we must count how many numbers $n \le L$ satisfy that condition when $T=100$ and $L=10^9$. If we write the set of admissible numbers as $\mathcal{H}(T,L)$, then the task is to determine $|\mathcal{H}(100,10^9)|$. The key observation is that only the primes up to $T$ matter: once those primes are fixed, every allowed composite factor is already built from them. Mathematical Approach Let $p_1 \lt p_2 \lt \cdots \lt p_r$ be the primes not exceeding $T$. The implementations count generalized Hamming numbers by exploring the exponent choices of these primes in increasing order. Prime Factors Are the Entire State Space By the fundamental theorem of arithmetic, every admissible number has a unique representation $$n=\prod_{i=1}^{r} p_i^{e_i}, \qquad e_i \ge 0,$$ and the condition $n \le L$ becomes $$\prod_{i=1}^{r} p_i^{e_i} \le L.$$ So the problem is not about testing arbitrary integers one by one. It is a structured counting problem over exponent vectors $(e_1,\dots,e_r)$. Composite numbers $\le T$ do not need their own branches, because their prime factors are already among the $p_i$....

Detailed mathematical approach

Problem Summary

A positive integer is a generalised Hamming number of type $n$ if none of its prime factors is larger than $n$. The Hamming numbers are those of type $5$, and $1$, which has no prime factor at all, is one of them: the statement lists the first Hamming numbers as $1, 2, 3, 4, 5, 6$, $8, 9, 10, 12, 15$, says that $1105$ of them do not exceed $10^8$, and asks how many generalised Hamming numbers of type $100$ do not exceed $10^9$. ^[1] In number theory they are called $n$-smooth, and $\Psi(x,y)$ counts the $y$-smooth integers up to $x$; ^[2] we are given $\Psi(10^8,5)=1105$ and must find $\Psi(10^9,100)$.

The whole solution in one sentence: every such number is a product of powers of the $25$ primes up to $100$ in exactly one way, so we choose the exponent of $2$, then that of $3$, and so on up to $97$, never letting the product exceed $10^9$, and count the complete choices.

Mathematical Approach

Throughout, $L$ is the limit and $p_1\lt p_2\lt\dots\lt p_r$ are the primes up to the type; for type $100$ these are the $r=25$ primes from $2$ to $97$.

1. Small cases

Up to $30$ there are $18$ Hamming numbers: $1, 2, 3, 4, 5, 6$, $8, 9, 10, 12, 15, 16$, $18, 20, 24, 25, 27, 30$. The other twelve, such as $7$ and $28=4\cdot7$, have a prime factor larger than $5$. Up to $10^0, 10^1,\dots,10^8$ there are $1, 9, 34, 86, 175$, $313, 507, 768, 1105$ Hamming numbers. ^[3] For type $100$ all numbers up to $100$ qualify, but then the share drops fast: $\Psi(10^3,100)=665$, $\Psi(10^4,100)=3716$, $\Psi(10^5,100)=17\,442$ and $\Psi(10^6,100)=72\,271$. The numbers we want are rare, so we build them instead of testing every integer up to $10^9$.

2. Only the primes up to 100 matter

By the fundamental theorem of arithmetic every integer greater than $1$ is a product of primes in exactly one way, up to the order of the factors, and $1$ is the empty product. ^[4] So the numbers we count correspond one to one to the vectors $(e_1,\dots,e_r)$ of non-negative integer exponents with

$$p_1^{e_1}p_2^{e_2}\cdots p_r^{e_r}\le L.$$

The primes come from the sieve of Eratosthenes: ^[5] every $i=2,\dots,10$ that is still unmarked is a prime and marks $i^2, i^2+i, \dots$ up to $100$; we may stop at $10$ because $11^2\gt100$. A composite $m\le100$ is marked by its smallest prime factor $q$, since $m=q\cdot(m/q)$ with $m/q\ge q$ gives $q^2\le m$, and a prime is never marked, as only multiples $ji$ with $j\ge i\ge2$ are. The $25$ unmarked numbers are the primes $2, 3, 5, \dots, 89, 97$.

3. Choosing the exponents one prime at a time

We fix the exponents from left to right. Let $F(i,m)$ be the number of integers $1\le u\le m$ whose prime factors all lie among $p_i,\dots,p_r$. Each such $u$ is $p_i^k v$ with a unique $k\ge0$ and a cofactor $v$ whose prime factors lie among $p_{i+1},\dots,p_r$, and since $v$ is an integer, $p_i^k v\le m$ holds exactly when $v\le\lfloor m/p_i^k\rfloor$. Sorting the numbers $u$ by $k$ gives

$$F(i,m)=\sum_{\substack{k\ge0\\ p_i^k\le m}}F\bigl(i+1,\lfloor m/p_i^k\rfloor\bigr),$$

with $F(r+1,m)=1$ for $m\ge1$, since only $u=1$ uses no prime. We want $F(1,L)$.

4. Worked example: the Hamming numbers up to 30

For type $5$ we have $p_1=2$, $p_2=3$ and $p_3=5$. Sorting the Hamming numbers up to $30$ by their power of two, $2^a=1, 2, 4, 8, 16$, leaves the bounds $\lfloor 30/2^a\rfloor=30, 15, 7, 3, 1$ for the odd part:

$a=0$: the odd parts $1, 3, 5, 9, 15, 25, 27$, so $F(2,30)=7$;
$a=1$: the odd parts $1, 3, 5, 9, 15$ give $2, 6, 10, 18, 30$, so $F(2,15)=5$;
$a=2$: the numbers $4, 12, 20$, so $F(2,7)=3$;
$a=3$: the numbers $8, 24$, so $F(2,3)=2$;
$a=4$: the number $16$, so $F(2,1)=1$.

So $F(1,30)=7+5+3+2+1=18$, and these are the numbers of Section 1. One level down, the powers $1, 3, 9, 27$ of three leave the bounds $30, 10, 3, 1$ for the powers $1, 5, 25$ of five, so $F(2,30)=3+2+1+1=7$.

5. The search tree

The programs walk through the tree of this recursion by depth-first search, following one branch to the bottom before trying the next. ^[6] Instead of the bound $m$ they carry the product $x$ of the prime powers chosen so far; the bound is $\lfloor L/x\rfloor$. A call at depth $i$ has fixed the exponents of $p_1,\dots,p_i$ and holds $x\le L$. At depth $r$ it counts one number; otherwise, with $p=p_{i+1}$, it calls depth $i+1$ with $x, xp, xp^2, \dots$ while the value does not exceed $L$. For type $5$ and $L=30$ the first call holds $x=1$ and calls depth $1$ with $1, 2, 4, 8, 16$, the branches of Section 4.

Claim. The calls at depth $i$ hold exactly the numbers up to $L$ whose prime factors lie among $p_1,\dots,p_i$, each of them once. So the calls at depth $r$ correspond one to one to the numbers we want.

Proof. By induction on $i$; depth $0$ holds only $x=1$. Every value $xp^k\le L$ passed on from depth $i$, where $p=p_{i+1}$, has its prime factors among $p_1,\dots,p_{i+1}$. Conversely, let $n\le L$ be such a number, $k$ the exponent of $p$ in $n$ and $x=n/p^k$. Then $x\le L$ is held at depth $i$, and because $xp^k=n\le L$, that call reaches the exponent $k$ and passes $n$ on. Since $n$ determines $k$ and $x$, and $x$ is held only once, $n$ arrives only once. $\square$

6. The stopping test

To decide whether $xp$ still fits, the programs compare $x$ with $\lfloor L/p\rfloor$ and never form a product above $L$. For positive integers the two tests agree:

$$xp\le L\iff x\le\lfloor L/p\rfloor,$$

because $x\le\lfloor L/p\rfloor$ gives $xp\le p\lfloor L/p\rfloor\le L$, and $xp\le L$ gives $x\le L/p$, so the integer $x$ is at most $\lfloor L/p\rfloor$. The first call of the example stops after $16$, as $16\gt\lfloor30/2\rfloor=15$. Every product the search forms is therefore at most $L$, and so is the count, so the unsigned 64-bit C++ search cannot overflow for any limit up to $2^{64}-1$.

7. The size of the tree

Put $p_0=1$. By the claim, depth $i$ has $\Psi(L,p_i)$ calls; for type $5$ and $L=30$ that makes $1+5+12+18=36$ calls. A counted number $n\gt1$ whose largest prime factor is $p_j$ is held at every depth from $j$ to $r$, and $n=1$ at every depth, so the total is

$$\sum_{i=0}^{r}\Psi(L,p_i)\le(r+1)\,\Psi(L,p_r).$$

Most calls only pass a finished number down: for type $100$ there are $721\,908$ calls for the $72\,271$ numbers up to $10^6$, and about $7.2$ calls per number at the real limit. A call whose next product $xp$ exceeds $L$ could count $1$ and return at once, since all later primes are larger; the programs do not take this shortcut. With it, every call that goes on would have at least two children, for the exponents $0$ and $1$, so a tree with $I$ such calls and $K$ leaves would have $I+K-1\ge2I$ edges, hence $I\lt K$: fewer than $2$ calls per counted number, about $1.9$ at the real limit.

8. Estimates and other methods

Taking logarithms turns the condition of Section 2 into $\sum_i e_i\ln p_i\le\ln L$. For type $5$ its real solutions with non-negative coordinates form a tetrahedron with perpendicular edges $\ln L/\ln2$, $\ln L/\ln3$ and $\ln L/\ln5$, whose volume is one sixth of their product, and the exponent vectors are its integer points. Rounding down the coordinates of a point of the tetrahedron keeps it inside, so the unit cubes with these points as lowest corners cover it, and the count is at least the volume $(\ln L)^3/(6\ln2\,\ln3\,\ln5)$: $1105$ against about $850$ for $L=10^8$. The refined formula

$$\frac{\bigl(\ln(L\sqrt{30})\bigr)^3}{6\ln2\,\ln3\,\ln5}$$

differs from the count by $O(\log L)$ and gives about $1107.8$ at $L=10^8$. ^[7] For $r$ primes the tetrahedron becomes an $r$-dimensional simplex, and for fixed primes the count grows like its volume $(\ln L)^r/(r!\,\ln p_1\cdots\ln p_r)$, ^[2] but for $25$ primes and $L=10^6$ this is about $10^{-9}$ instead of $72\,271$, so we need an exact count.

Two other exact methods are worth knowing. A direct scan that divides every integer up to $10^9$ by the $25$ primes took about 3 seconds in C, a hundred times as long as the search. Dijkstra's method outputs $1$ and merges the sequences $2H$, $3H$ and $5H$ of the Hamming numbers $H$, emitting a value such as $6=2\cdot3=3\cdot2$ once; ^[7] it lists the numbers in increasing order, which counting does not need.

How the Code Works

The C++ program ^[8] keeps the limit $10^9$, the type $100$ and a checkpoint switch in Options. parse_arguments accepts --limit=, --type= and --skip-checkpoints; its digit readers parse_u64_after_prefix and parse_int_after_prefix do not check for overflow, and a type close to $2^{31}$ would also overflow the int arithmetic of the sieve. An unknown argument, a limit below $1$ or a type below $2$ gives exit code $1$. run_checkpoints requires count_hamming(100000000ULL, 5) to return the $1105$ of the statement, or the program stops with exit code $2$; then main prints the count.

primes_up_to is the sieve of Section 2, with one byte per number in a std::vector<std::uint8_t>. count_hamming runs the search of Sections 5 and 6 as a recursive lambda dfs, which receives itself as self, a common way to let a lambda declared with auto call itself. Its idx counts from $0$, so primes[idx] is $p_{i+1}$ at depth $i$, and value is $x$. When idx reaches the number of primes, the call increments the 64-bit count; otherwise its loop calls the next depth, tests x > limit / p and multiplies by p only if that test fails.

The Python and Java programs fix the limit and the type and have no options or checkpoint. Python sieves a list of booleans, and its nested function dfs updates count through nonlocal, with the test x > limit // p. Java turns its list of primes into an int[] and recurses in a static method dfs. C++ uses unsigned and Java signed 64-bit integers, which suffice because at the problem's size no value exceeds $10^9$, and Python's integers have no fixed size. ^[8] For type $100$ the recursion is $26$ calls deep, in general $r+1$.

Complexity and Verification

The sieve takes $O(T\log\log T)$ steps for the type $T=100$. ^[5] By Section 7 the search makes at most $(r+1)\,\Psi(L,T)$ calls, and each loop step makes one child call, so there are fewer loop steps than calls. With arithmetic counted as constant time, the total time is $O(T\log\log T+r\,\Psi(L,T))$, which depends on $L$ only through the count, and the memory is $O(T+r)$: the sieve, the primes and a stack of $r+1$ calls. On an Apple M4 Max laptop, measured with the system's time command, a full run takes about $0.03$ s in C++ compiled with -O2, including the checkpoint, $0.05$ s in Java and $1.0$ s in Python.

The published code checks only $\Psi(10^8,5)=1105$, in C++. For this article the C++ program was run at $340$ pairs of a limit up to $200\,000$ and a type from $2$ to $150$, and copies of the Python and Java programs with the two constants changed at $7$ pairs; all agreed with a count of the integers whose largest prime factor does not exceed the type. A Python check on $153$ pairs confirmed the call counts per depth of Section 7. At the real parameters the three programs, the search with the shortcut and an independent C program that divides every integer up to $10^9$ by the primes up to $97$ all give the same count. These checks are not part of the published solutions.

Footnotes and References

The references below supply the background facts; the claims about the search tree and its size are derived in this article.

Project Euler 204 — Generalised Hamming Numbers. The official statement defines Hamming numbers as positive numbers with no prime factor larger than 5, lists the first eleven, states that 1105 of them do not exceed 10⁸, and asks how many numbers with no prime factor larger than 100 do not exceed 10⁹.
Wikipedia — Smooth number. A positive integer is B-smooth if none of its prime factors is greater than B. The de Bruijn function Ψ(x, y) counts the y-smooth integers up to x, and for a fixed small B, Ψ(x, B) ~ (1/π(B)!) ∏ log x / log p over the primes p ≤ B.
OEIS A106598. Number of positive integers ≤ 10ⁿ divisible by no prime exceeding 5: 1, 9, 34, 86, 175, 313, 507, 768, 1105, … for n = 0, 1, 2, …
Wikipedia — Fundamental theorem of arithmetic. Every integer greater than 1 is a prime or a product of primes, unique up to the order of the factors; by convention the empty product is 1.
Wikipedia — Sieve of Eratosthenes. Marking the multiples of a prime p may start at p², the sieve may stop once p² exceeds n, and it takes O(n log log n) operations.
Wikipedia — Depth-first search. Depth-first search explores as far as possible along each branch before backtracking.
Wikipedia — Regular number. Hamming numbers, the 5-smooth numbers, are named after Richard Hamming, who proposed generating them in ascending order; Edsger Dijkstra popularised algorithms for this, such as merging 2H, 3H and 5H. Their number up to N is about the volume log₂N · log₃N · log₅N / 6 of a tetrahedron, and more precisely (ln(N√30))³/(6 ln 2 ln 3 ln 5) + O(log N).
C++ — Euler204.cpp. The linked immutable revision contains parse_arguments, primes_up_to, count_hamming and run_checkpoints; it is the implementation record for this article. Python — numeric types: integers have unlimited precision.

Problem 204 source code

C++

#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using u64 = std::uint64_t;

struct Options {
    u64 limit = 1000000000ULL;
    int type = 100;
    bool run_checkpoints = true;
};

bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10ULL + static_cast<u64>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);

        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_u64_after_prefix(arg, "--limit=", options.limit) ||
            parse_int_after_prefix(arg, "--type=", options.type)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.limit >= 1 && options.type >= 2;
}

std::vector<int> primes_up_to(const int n) {
    std::vector<std::uint8_t> is_prime(static_cast<std::size_t>(n + 1), 1);
    is_prime[0] = 0;
    is_prime[1] = 0;

    for (int i = 2; i * i <= n; ++i) {
        if (!is_prime[static_cast<std::size_t>(i)]) {
            continue;
        }
        for (int j = i * i; j <= n; j += i) {
            is_prime[static_cast<std::size_t>(j)] = 0;
        }
    }

    std::vector<int> primes;
    for (int i = 2; i <= n; ++i) {
        if (is_prime[static_cast<std::size_t>(i)]) {
            primes.push_back(i);
        }
    }
    return primes;
}

u64 count_hamming(const u64 limit, const int type) {
    const std::vector<int> primes = primes_up_to(type);

    u64 count = 0;

    const auto dfs = [&](auto&& self, std::size_t idx, u64 value) -> void {
        if (idx == primes.size()) {
            ++count;
            return;
        }

        const u64 p = static_cast<u64>(primes[idx]);
        for (u64 x = value;;) {
            self(self, idx + 1, x);
            if (x > limit / p) {
                break;
            }
            x *= p;
        }
    };

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

bool run_checkpoints() {
    if (count_hamming(100000000ULL, 5) != 1105ULL) {
        std::cerr << "Checkpoint failed for type=5, limit=1e8" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << count_hamming(options.limit, options.type) << '\n';
    return 0;
}

Python

# Problem 204: Generalised Hamming Numbers
# Count 100-smooth numbers (Hamming) up to 10^9.

def solve():
    limit = 10**9
    type_n = 100
    # Sieve primes up to type_n
    sieve = [True] * (type_n + 1)
    sieve[0] = sieve[1] = False
    for i in range(2, int(type_n**0.5)+1):
        if sieve[i]:
            for j in range(i*i, type_n+1, i):
                sieve[j] = False
    primes = [i for i in range(2, type_n+1) if sieve[i]]

    count = 0
    def dfs(idx, val):
        nonlocal count
        if idx == len(primes):
            count += 1
            return
        p = primes[idx]
        x = val
        while True:
            dfs(idx + 1, x)
            if x > limit // p:
                break
            x *= p
    dfs(0, 1)
    print(count)

solve()

Java

public class Euler204 {
    static long limit = 1000000000L;
    static long count = 0;

    public static void main(String[] args) {
        int type = 100;
        boolean[] sieve = new boolean[type + 1];
        java.util.Arrays.fill(sieve, true);
        sieve[0] = sieve[1] = false;
        for (int i = 2; i * i <= type; i++)
            if (sieve[i])
                for (int j = i * i; j <= type; j += i)
                    sieve[j] = false;
        java.util.List<Integer> primes = new java.util.ArrayList<>();
        for (int i = 2; i <= type; i++)
            if (sieve[i])
                primes.add(i);
        int[] p = primes.stream().mapToInt(Integer::intValue).toArray();
        dfs(p, 0, 1);
        System.out.println(count);
    }

    static void dfs(int[] primes, int idx, long val) {
        if (idx == primes.length) {
            count++;
            return;
        }
        long p = primes[idx];
        for (long x = val;;) {
            dfs(primes, idx + 1, x);
            if (x > limit / p)
                break;
            x *= p;
        }
    }
}