Problem 187: Semiprimes

Project Euler Problem 187 Solution

EulerSolve provides an optimized solution for Project Euler Problem 187, Semiprimes, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Project Euler 187 asks for the number of composite integers $n \lt 10^8$ that have exactly two prime factors when multiplicity is counted. In other words, we must count the semiprimes below the limit, including prime squares such as $4=2^2$, $9=3^2$, and $25=5^2$. Instead of factoring every integer below the limit, the solution counts prime pairs directly. If a number is semiprime, then it can be written as $n=pq$ with primes $p \le q$. Once that order is enforced, every valid semiprime corresponds to exactly one pair $(p,q)$. Mathematical Approach Let $L$ denote the upper bound. The target quantity is $$S(L)=\#\{n \lt L : n=pq,\ p \le q,\ p,q\text{ prime}\}.$$ The key idea is to choose the smaller prime first and then count how many larger primes can be paired with it. Ordered prime pairs give a one-to-one model Every semiprime has a factorization $n=pq$ with primes $p \le q$. Because the smaller factor is placed first, the representation is unique. So the whole problem becomes a counting problem on ordered prime pairs satisfying $$pq \lt L,\qquad p \le q.$$ This removes the need for individual factorization. We never inspect a candidate integer and ask whether it is semiprime; we build exactly the semiprimes we want by counting admissible pairs....

Detailed mathematical approach

Problem Summary

Write $\Omega(n)$ for the number of prime factors of $n$ counted with multiplicity. The examples of the statement are $15=3\cdot5$ and $9=3\cdot3$, with $\Omega=2$, and $12=2\cdot2\cdot3$, with $\Omega(12)=3$. We must count the composite integers $n\lt10^8$ that have exactly two prime factors, not necessarily distinct, that is, the $n\lt10^8$ with $\Omega(n)=2$. These numbers are called semiprimes. ^[2] The statement lists the ten of them below $30$: $4, 6, 9, 10, 14$, $15, 21, 22, 25, 26$. ^[1]

The whole solution in one sentence: every such number is $pq$ for exactly one pair of primes $p\le q$, so we list all primes up to $(10^8-1)/2$ with the sieve of Eratosthenes and, for each prime $p$ with $p^2\lt10^8$, count the primes $q$ from $p$ up to $(10^8-1)/p$ by a binary search in that sorted list.

Mathematical Approach

Write $L$ for the limit, a positive integer, and $S(L)$ for the number of positive integers $n\lt L$ with $\Omega(n)=2$. The problem asks for $S(10^8)$, and the statement gives $S(30)=10$.

1. Small cases

Sorting the numbers from $2$ to $29$ by $\Omega$ gives three groups:

$\Omega=1$, the primes: $2, 3, 5, 7, 11$, $13, 17, 19, 23, 29$;
$\Omega=2$: $4, 6, 9, 10, 14$, $15, 21, 22, 25, 26$;
$\Omega\ge3$: $8, 12, 16, 18$, $20, 24, 27, 28$.

The number $1$ has no prime factor at all. The middle group is the list of the statement, so $S(30)=10$. Counting further gives $S(10)=3$, $S(100)=34$, $S(10^3)=299$, $S(10^4)=2625$, $S(10^5)=23378$, $S(10^6)=210035$ and $S(10^7)=1904324$. Rather than factor every number below $10^8$, we build the semiprimes from the primes.

2. One semiprime, one pair of primes

Claim. With $\mathbb P$ the set of primes, the rule $(p,q)\mapsto pq$ matches the pairs of primes with $p\le q$ and $pq\lt L$ one to one with the integers $n\lt L$ with $\Omega(n)=2$. Hence

$$\begin{aligned}S(L)=\#\{(p,q)\in\mathbb P^2:\ &p\le q,\\&pq\lt L\}.\end{aligned}$$

Proof. By the fundamental theorem of arithmetic every integer $n\ge2$ is a product of primes in exactly one way, up to the order of the factors. ^[3] For primes $p$ and $q$ the factorisation of $pq$ is $p\cdot q$, so $\Omega(pq)=2$. Conversely, if $\Omega(n)=2$, the factorisation of $n$ consists of exactly two primes, and listing the smaller one first gives a pair with $p\le q$ and $pq=n$. If two such pairs $(p,q)$ and $(p',q')$ have $pq=p'q'$, they are two factorisations of the same number, so they consist of the same two primes, and since both pairs are in nondecreasing order, $(p,q)=(p',q')$. $\square$

The order $p\le q$ makes the count exact: $6$ comes only from $(2,3)$, not also from $(3,2)$, and the square $25$ from the single pair $(5,5)$.

3. Where the pairs live

If $p\le q$ and $pq\lt L$, then $p^2\le pq\lt L$, so the smaller prime satisfies $p^2\lt L$. Since $pq$ and $L$ are integers, $pq\lt L$ means $pq\le L-1$, and because $p\ge2$,

$$q\le\Bigl\lfloor\frac{L-1}{p}\Bigr\rfloor\le\Bigl\lfloor\frac{L-1}{2}\Bigr\rfloor=M.$$

A list of the primes up to $M$ therefore contains every prime that occurs in a pair. For $L=10^8$ we get $M=49\,999\,999$, and the smaller prime runs over the $1229$ primes below $10^4$, up to $9973$; ^[5] the next prime, $10007$, already has $10007^2\gt10^8$.

4. Counting the partners of one prime

Let $p_1=2$, $p_2=3$, $p_3=5,\dots$ be the primes in increasing order, and let $\pi(x)$, the prime-counting function, be the number of primes $\le x$, so that $\pi(p_k)=k$. ^[5] For the smaller prime $p=p_k$, the partners $q$ are the primes with $p_k\le q\le\lfloor(L-1)/p_k\rfloor$. There are $\pi(\lfloor(L-1)/p_k\rfloor)$ primes up to the upper end, and the $k-1$ primes $p_1,\dots,p_{k-1}$ below $p_k$ must be left out. Adding up over all possible smaller primes gives

$$\begin{aligned}S(L)=\sum_{k:\ p_k^2\lt L}\Bigl(&\pi\Bigl(\Bigl\lfloor\frac{L-1}{p_k}\Bigr\rfloor\Bigr)\\&{}-k+1\Bigr),\end{aligned}$$

which is the known formula for the number of semiprimes up to $L-1$. ^[2] Every term is at least $1$, because $p_k^2\lt L$ makes $p_k$ a partner of itself.

5. Worked example: $L=30$

Here $M=\lfloor29/2\rfloor=14$, and the primes up to $14$ are $2,3,5,7,11,13$. The smaller prime must satisfy $p^2\lt30$, so $p\in\{2,3,5\}$; for $p=7$ already $49\ge30$.

$p=2$: $q\le\lfloor29/2\rfloor=14$, so $q\in\{2,3,5,7,11,13\}$, giving $4,6,10,14,22,26$; the term is $\pi(14)-1+1=6$.
$p=3$: $q\le\lfloor29/3\rfloor=9$, so $q\in\{3,5,7\}$, giving $9,15,21$; the term is $\pi(9)-2+1=3$.
$p=5$: $q\le\lfloor29/5\rfloor=5$, so $q=5$, giving $25$; the term is $\pi(5)-3+1=1$.

The total is $6+3+1=10$, and the ten products are exactly the list of the statement. For $L=100$ the bounds for $p=2,3,5,7$ are $49$, $33$, $19$ and $14$, and the terms $15+10+6+3=34$ agree with Section 1.

6. Listing the primes: the sieve of Eratosthenes

The sieve goes through the numbers $2,3,\dots,M$ in increasing order. ^[4] Whenever the current number $p$ is unmarked, it is a prime, and we mark its multiples $p^2$, $p^2+p$, $p^2+2p,\dots$ up to $M$. We stop as soon as $p^2\gt M$. For $M=14$, the prime $2$ marks $4,6,8,10,12,14$, the prime $3$ marks $9$ and $12$, and at $5$ already $25\gt14$; the unmarked numbers $2,3,5,7,11,13$ are the primes up to $14$.

Why this is correct: a marked number is a multiple $jp$ with $j\ge p\ge2$, so no prime is marked. A composite $n\le M$ with smallest prime factor $r$ is $n=r\cdot(n/r)$ with $n/r\ge r$, so $r^2\le n\le M$; the prime $r$ is therefore used, and $n$ is one of the numbers $r^2,r^2+r,\dots$ that it marks, so starting at $r^2$ loses nothing. For $M=49\,999\,999$ the marking primes are the $908$ primes up to $\lfloor\sqrt M\rfloor=7071$, and the sieve leaves $\pi(M)=3\,001\,134$ primes.

7. What the sieve costs

The prime $p$ makes $\lfloor(M-p^2)/p\rfloor+1$ marking steps, fewer than $M/p$. Summing over the primes up to $\sqrt M$ and using Mertens' second theorem, ^[6]

$$\sum_{p\le x}\frac1p=\ln\ln x+B+O\Bigl(\frac{1}{\ln x}\Bigr),$$

where $B=0.2614972\ldots$ is the Meissel–Mertens constant, shows that the sieve takes $O(M\log\log M)$ steps. ^[4] For $M=49\,999\,999$ it makes $119\,266\,369$ marking steps, about $2.39$ per entry, and the reciprocals of the $908$ marking primes add up to about $2.4452$, close to $\ln\ln7071+B\approx2.4435$.

8. Counting the partners by binary search

The sieve produces the primes in increasing order, $P[0]\lt P[1]\lt\dots$, with $P[i]=p_{i+1}$. In a sorted list of $n$ entries, numbered from $0$, let $j$ be the index of the first entry $\gt x$, or $j=n$ if there is none; then exactly $j$ entries are $\le x$. Binary search finds $j$ by comparing $x$ with the middle entry of the remaining range and keeping the half that can contain the position; it needs at most $\lfloor\log_2n\rfloor+1$ comparisons, ^[7] that is $22$ for the $3\,001\,134$ primes up to $M$. For the smaller prime $p=P[i]$ we search with $x=\lfloor(L-1)/p\rfloor$: the partners occupy the indices $i$ to $j-1$, so their number is $j-i$, the term of Section 4 with $k=i+1$. For $L=30$ the list is $2,3,5,7,11,13$; for $p=3$ we have $i=1$, $x=9$ and $j=4$, the index of $11$, so the partners are the $3$ primes $3,5,7$, and for $p=2$ no entry exceeds $14$, so $j=6$. At the real limit the last smaller prime, $9973$, has the bound $\lfloor(10^8-1)/9973\rfloor=10027$ and the three partners $9973$, $10007$ and $10009$.

All $1229$ searches make at most $1229\cdot22=27038$ comparisons, nothing next to the sieve. Because the bounds shrink as $p$ grows, a pointer moving only down the list would do as well. For far larger limits, sieving up to $L/2$ would give way to the combinatorial method of Meissel and Lehmer for the values $\pi(x)$. ^[5]

How the Code Works

The C++ program ^[8] keeps the limit $L=10^8$ and a checkpoint switch in Options. parse_arguments accepts --skip-checkpoints and --limit= with decimal digits, which parse_int_after_prefix reads into an int without an overflow check, so a limit of $2^{31}$ or more overflows the int undetected; the limit must stay below $2^{31}$. Any other argument, or a limit below $4$, gives exit code $1$. Unless switched off, run_checkpoints requires solve(30) to return $10$ and solve(10000) to equal brute_small(10000), which factors every $n$ from $4$ to $9999$ by trial division, stopping at a third prime factor; a failure gives exit code $2$. Then main prints solve(options.limit).

sieve_primes is the sieve of Section 6: a std::vector<std::uint8_t> named is_prime holds one byte per number, $1$ for unmarked and $0$ for marked, with the numbers $0$ and $1$ marked from the start. Every p that is still unmarked and whose square, formed in 64 bits, does not exceed the argument marks p * p, p * p + p and so on, and the unmarked indices go into a std::vector<int>. solve sieves up to max_q = (limit - 1) / 2, the $M$ of Section 3. For each index i it takes p = primes[i] as a 64-bit integer, stops when p * p >= limit, computes q_max = (limit - 1) / p and calls std::upper_bound from primes.begin() + i to the end, which returns the first prime greater than q_max, or the end of the list. ^[9] Its distance from position i is the term of Section 4, added to a long long total.

The Python and Java programs fix the limit at $10^8$, without options or checkpoints. Python sieves in a bytearray with an inner loop over range(i*i, max_q+1, i) for every unmarked i up to int(max_q**0.5), which is $7071$, and calls bisect.bisect_right(primes, qmax) on the whole list; ^[9] subtracting i removes the primes below $p$. Java sieves in a boolean[] and stores the primes in an int[] of length maxQ, room for $49\,999\,999$ entries of which $3\,001\,134$ are used. Its own binary search starts with lo = i, hi = pc - 1 and ans = i - 1, moves ans to the last index whose prime is at most qmax and adds ans - i + 1. The count is below $L$, so the 64-bit totals of C++ and Java cannot overflow.

Complexity and Verification

With $M=\lfloor(L-1)/2\rfloor$, the sieve takes $O(M\log\log M)$ time by Section 7, collecting the primes takes $O(M)$, and the $\pi(\sqrt{L-1})$ binary searches take $O(\log M)$ each, so the sieve dominates. The memory is one byte per number up to $M$, $50$ MB, plus the primes: in C++ about $12$ MB of 4-byte values in a std::vector that may hold some spare capacity, while the byte array is freed when sieve_primes returns; in Python a list of integer objects; in Java the whole $200$ MB array. Measured with the system's time command on a laptop with an Apple M4 Max chip, using Apple clang with -O2, OpenJDK 17 and CPython 3.11, a full run takes about $0.1$ s in C++, $0.2$ s in Java and $2.6$ s in Python, with a peak resident memory of about $81$, $297$ and $186$ MB.

The C++ checkpoints compare solve(30) with the $10$ of the statement and solve(10000) with the brute-force count $2625$. For this article the C++ program was run with --limit=L for every $L$ from $4$ to $3000$, agreeing each time with a count by trial division, and a Python version of the formula of Section 4 agreed with a sieve for $\Omega(n)$ for every $L$ from $4$ to $20000$. Copies of the Python and Java programs with only the fixed limit changed gave the trial-division counts at $16$ limits from $4$ to $20000$. An independent C program computed $\Omega(n)$ for every $n\lt10^8$ with a segmented sieve and counted the $n$ with $\Omega(n)=2$, without forming any prime pairs; it gave the same result as the three published programs. These checks are not part of the published solutions.

Footnotes and References

The references below supply the background facts used above; the bounds of Section 3, the worked examples and the operation counts are derived in this article.

Project Euler 187 — Semiprimes. The official statement lists the ten composites below thirty with precisely two, not necessarily distinct, prime factors and asks how many composite integers n < 10⁸ have precisely two such factors.
Wikipedia — Semiprime. A semiprime is a product of exactly two primes, which may be equal. The article gives the number of semiprimes up to n as the sum, over k from 1 to π(√n), of π(n/pₖ) − k + 1, where pₖ, the k-th prime, is the smaller factor.
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.
Wikipedia — Sieve of Eratosthenes. Marking the multiples of p may start at p², and the process may stop once p² exceeds n. The algorithm needs O(n log log n) operations, because the sum of the reciprocals of the primes grows like log log n.
Wikipedia — Prime-counting function. π(x) counts the primes up to x, and the table gives π(10⁴) = 1229. A sieve finds π(x) if x is not too large; Meissel described a practical combinatorial way of evaluating π(x), which Lehmer extended in 1959.
Wikipedia — Mertens' theorems. The second theorem: the sum of 1/p over the primes p ≤ n is log log n + M + O(1/log n), where M ≈ 0.2614972128, written B in this article, is the Meissel–Mertens constant.
Wikipedia — Binary search. Binary search compares the target with the middle element of a sorted array and discards the half that cannot contain it; in the worst case it makes ⌊log₂ n⌋ + 1 iterations.
C++ — Euler187.cpp. The linked immutable revision contains parse_arguments, sieve_primes, solve, brute_small and run_checkpoints. It is the implementation record for this article.
cppreference — std::upper_bound returns the first element of a sorted range that is ordered after the given value, with at most log₂ N + O(1) comparisons. In Python's bisect module, bisect_right returns an insertion point with every element to its left at most x.

Problem 187 source code

C++

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

namespace {

using i64 = long long;

struct Options {
    int limit = 100000000;
    bool run_checkpoints = 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_int_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.limit >= 4;
}

std::vector<int> sieve_primes(const int limit) {
    std::vector<std::uint8_t> is_prime(static_cast<std::size_t>(limit + 1), 1U);
    if (limit >= 0) {
        is_prime[0] = 0U;
    }
    if (limit >= 1) {
        is_prime[1] = 0U;
    }
    for (int p = 2; static_cast<i64>(p) * p <= limit; ++p) {
        if (is_prime[static_cast<std::size_t>(p)] == 0U) {
            continue;
        }
        for (int q = p * p; q <= limit; q += p) {
            is_prime[static_cast<std::size_t>(q)] = 0U;
        }
    }
    std::vector<int> primes;
    for (int p = 2; p <= limit; ++p) {
        if (is_prime[static_cast<std::size_t>(p)] != 0U) {
            primes.push_back(p);
        }
    }
    return primes;
}

i64 solve(const int limit) {
    const int max_q = (limit - 1) / 2;
    const std::vector<int> primes = sieve_primes(max_q);

    i64 count = 0;
    const int total = static_cast<int>(primes.size());
    for (int i = 0; i < total; ++i) {
        const i64 p = primes[static_cast<std::size_t>(i)];
        if (p * p >= limit) {
            break;
        }
        const i64 q_max = (limit - 1) / p;
        const auto it = std::upper_bound(primes.begin() + i, primes.end(), static_cast<int>(q_max));
        count += static_cast<i64>(it - (primes.begin() + i));
    }
    return count;
}

i64 brute_small(const int limit) {
    int count = 0;
    for (int n = 4; n < limit; ++n) {
        int x = n;
        int omega = 0;
        for (int p = 2; p * p <= x; ++p) {
            while (x % p == 0) {
                x /= p;
                ++omega;
                if (omega > 2) {
                    break;
                }
            }
            if (omega > 2) {
                break;
            }
        }
        if (x > 1) {
            ++omega;
        }
        if (omega == 2) {
            ++count;
        }
    }
    return count;
}

bool run_checkpoints() {
    if (solve(30) != 10) {
        std::cerr << "Checkpoint failed for limit 30 sample" << '\n';
        return false;
    }
    if (solve(10000) != brute_small(10000)) {
        std::cerr << "Checkpoint failed for brute cross-check at limit 10000" << '\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 << solve(options.limit) << '\n';
    return 0;
}

Python

# Problem 187: Semiprimes
# Count composites below 10^8 that are the product of exactly two primes.

def solve():
    limit = 10**8
    max_q = (limit-1)//2
    # Sieve of Eratosthenes
    sieve = bytearray(b'\x01') * (max_q+1)
    sieve[0] = 0
    sieve[1] = 0
    for i in range(2, int(max_q**0.5)+1):
        if sieve[i]:
            for j in range(i*i, max_q+1, i):
                sieve[j] = 0
    primes = [i for i in range(2, max_q+1) if sieve[i]]

    import bisect
    count = 0
    for i, p in enumerate(primes):
        if p*p >= limit: break
        qmax = (limit-1)//p
        j = bisect.bisect_right(primes, qmax)
        count += j - i
    print(count)

solve()

Java

public class Euler187 {
    public static void main(String[] args) {
        int limit = 100000000, maxQ = (limit - 1) / 2;
        boolean[] sieve = new boolean[maxQ + 1];
        java.util.Arrays.fill(sieve, true);
        sieve[0] = false;
        sieve[1] = false;
        for (int i = 2; (long) i * i <= maxQ; i++)
            if (sieve[i])
                for (int j = i * i; j <= maxQ; j += i)
                    sieve[j] = false;
        int[] primes = new int[maxQ];
        int pc = 0;
        for (int i = 2; i <= maxQ; i++)
            if (sieve[i])
                primes[pc++] = i;
        long count = 0;
        for (int i = 0; i < pc; i++) {
            long p = primes[i];
            if (p * p >= limit)
                break;
            long qmax = (limit - 1) / p;
            int lo = i, hi = pc - 1, ans = i - 1;
            while (lo <= hi) {
                int mid = (lo + hi) / 2;
                if (primes[mid] <= qmax) {
                    ans = mid;
                    lo = mid + 1;
                } else
                    hi = mid - 1;
            }
            count += ans - i + 1;
        }
        System.out.println(count);
    }
}