Project Euler 7: The 10001st Prime with a Sieve of Eratosthenes That Doubles in Size

This is the seventh problem in the Project Euler series. It asks us to count primes in order, so the trial division from problem 3 already gets you an answer. This post goes a step further: it derives the sieve of Eratosthenes and answers the question every sieve raises, “how far do I sieve?”, not with an unexplained constant but by doubling the sieve until it is big enough. The sieve comes back in later problems such as problem 10.

Spoiler warning: this post contains a full solution and the answer. If you want to try it yourself first, start with the original problem. Posts in this series never carry ads, because the problem text is licensed under CC BY-NC-SA 4.0, which is non-commercial.

Problem

Listing the primes from smallest up, the first six are 2, 3, 5, 7, 11 and 13, so the sixth prime is 13. Counting the same way, the task is to find the 10001st prime.

There is no input; n is fixed at 10001. The output is a single integer. In this post pₙ denotes the n-th prime (p₁ = 2, p₆ = 13).

Background

Primes and composites. An integer greater than 1 whose only divisors are 1 and itself is prime; every other integer greater than 1 is composite. 1 is neither, which is why the count starts at p₁ = 2.

Why dividing up to √m is enough. If m = a × b is composite with 1 < a ≤ b, then a × a ≤ a × b = m, so a ≤ √m. Every composite therefore has a divisor no larger than √m. Conversely, an m ≥ 2 that no number from 2 to √m divides is prime. In code we write d ≤ m / d (integer division) instead of d × d ≤ m, which says the same thing without a multiplication that could overflow: for positive integers, d ≤ ⌊m / d⌋ and d × d ≤ m are exactly the same condition.

The smallest divisor is prime. Let q be the smallest divisor greater than 1 of some m ≥ 2. Then q is prime: if it were composite it would have a smaller divisor greater than 1, which would also divide m, contradicting the choice of q. Both proofs below use this.

There is no last prime (Euclid). Suppose p₁, …, pₖ were all the primes. The smallest divisor greater than 1 of N = p₁ × … × pₖ + 1 is prime, but dividing N by any pᵢ leaves remainder 1, so that prime is not on the list. This contradiction shows there are infinitely many primes. It is the reason “keep growing the sieve and you will eventually collect n primes” holds, that is, why the algorithm below terminates.

The sieve of Eratosthenes. Write out the numbers 0 to L. Going up from 2, take each number that has not been crossed out, record it as prime, and cross out its multiples. When you are done, the numbers ≥ 2 left uncrossed are exactly the primes up to L. The derivation builds this step by step.

Brute force

Step through the numbers from 2 upwards, test each one by trial division, and stop at the n-th prime.

def is_prime(m):
    if m < 2:
        return False
    d = 2
    while d <= m // d:
        if m % d == 0:
            return False
        d += 1
    return True

count = 0
m = 1
while count < 10001:
    m += 1
    if is_prime(m):
        count += 1
# m is now the 10001st prime

This is correct and needs no bound in advance. Each m costs at most √m divisions, so testing every number up to pₙ costs at most about pₙ × √pₙ, that is, O(pₙ√pₙ). Composites stop early at a small divisor, but the n primes each need the full run up to √m, and they make up most of the cost. For n = 10001 this is still quick enough, but the factor √pₙ grows fast with n. Instead of asking each number “are you prime?”, we can start from the primes and cross out their multiples all at once, which removes almost all of that factor. That is the sieve.

Derivation

Step 1: check the problem’s example by hand. Counting primes from 2 gives 2 (1st), 3 (2nd), 5 (3rd), 7 (4th), 11 (5th), 13 (6th). The numbers in between, 4, 6, 8, 9, 10 and 12, are each divisible by 2 or 3. p₆ = 13 matches the problem’s example.

Step 2: cross out instead of asking. Trial division asks, for every m, “does some smaller number divide m?”. Turn it around: each time we find a prime p, mark all multiples of p as composite. Then, when we reach a number, the absence of a mark alone tells us it is prime. With L = 16:

i = 2: unmarked → prime. Cross out 4, 6, 8, 10, 12, 14, 16
i = 3: unmarked → prime. Cross out 9, 12, 15 (6 is already crossed out)
i = 4: crossed out → skip
i = 5: unmarked → prime. 5 × 5 = 25 > 16, nothing to cross out
i = 6 … 16: only 7, 11, 13 are unmarked → prime
Result: 2, 3, 5, 7, 11, 13

Six primes, and the sixth is 13, as in step 1.

Step 3: start crossing out at i × i, and skip it when i × i > L. Among the multiples k × i of a prime i, those with k < i are multiples of a prime factor of k, which is smaller than i, so they were crossed out when that prime was handled. Starting at i × i therefore misses nothing. Not crossing out 6 again at i = 3 above is an example. For the same reason, when i × i > L there is nothing left to cross out, and such an i is only checked and added to the list.

Step 4: a sieve needs a bound, and we don’t know how big pₙ is. A sieve gives “all primes up to L”, while the problem asks for “the n-th prime”. The common move is to pick some large L, but without a reason that L is large enough there is no guarantee the answer is right. This post does not fix L in advance:

• Sieve with L = 16.
• If that gives at least n primes, return the n-th of them.
• Otherwise double L and sieve again from scratch.

Why the n-th entry is the answer once there are at least n. The sieve returns the primes up to L all of them, in increasing order. So the n-th entry of the list is the n-th smallest prime, pₙ. That the list only covers numbers up to L does not matter: every prime smaller than pₙ is also smaller than pₙ ≤ L, so they are all on the list.

Why it terminates. There are infinitely many primes, so doubling L again and again eventually makes L ≥ pₙ, and at that point the list holds at least n primes.

Trace it on small cases. For n = 6, the first sieve at L = 16 already collects six primes, so it returns 13 at once. For n = 7, L = 16 has only six, which is not enough; sieving again at L = 32 gives eleven primes, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, and it returns the seventh, 17.

Step 5: doubling wastes little. Suppose the loop does not stop at the first L = 16, and let L be the size of the last sieve. The previous sieve at L / 2 held fewer than n primes, so pₙ > L / 2, and the last sieve succeeded, so pₙ ≤ L.

L / 2 < pₙ ≤ L,   so   L < 2pₙ

The last sieve is less than twice the size actually needed. Adding all the earlier sieves, 16 + 32 + … + L < 2L < 4pₙ, so compared with sieving once with an exact bound known in advance, this costs only a constant factor more. That is why it works without any proven bound on pₙ.

Final algorithm

1. primes_up_to(L): make a marker array of size L + 1. For i from 2 to L, if i is unmarked, append it to the list of primes, and if i ≤ L / i, cross out i × i, i × i + i, ….
2. nth_prime(n): starting from L = 16, repeat step 1, doubling L, until the list has at least n entries; return the n-th entry (index n − 1, counting from 0).
3. Before printing the answer, check: p₁ = 2 and p₆ = 13, and the first 100 primes found by trial division must equal nth_prime(1) through nth_prime(100). If anything differs, exit with a nonzero status.

Correctness proof

• Primes are never crossed out. Every crossed-out number has the form i × i + k × i = i × (i + k) with i ≥ 2 and k ≥ 0, a product of two numbers greater than 1, hence composite.
• Every composite up to L is crossed out. The smallest divisor q of a composite m ≤ L is prime, and by the background section q ≤ √m ≤ √L. When the loop reaches i = q, q is unmarked (it is prime) and q ≤ L / q, so it crosses out. Since m = q × (m / q) with m / q ≥ q, m lies in the range crossed out starting from q × q. And q < m, so this happens before the loop reaches m.
• So when the loop reaches i, “unmarked” is equivalent to “prime”, and since i runs through 2 to L in order, the list holds every prime up to L in increasing order.
• As shown in step 4, once the list has at least n entries its n-th entry is pₙ, and since there are infinitely many primes the loop terminates.

Complexity

Split the work of one sieve of size L into three parts.

• Building the sieve (crossing out). Each prime p ≤ √L crosses out at most L / p numbers. Even summing L / i over all i ≤ √L gives L × (1 + 1/2 + …) = O(L log L), an easy upper bound. Using the known result that the sum of the reciprocals of the primes grows only like log log (Mertens’ theorem), it is more precisely O(L log log L).
• Scanning (building the list). The outer loop checks each marker once for i = 2 to L: O(L).
• Selecting. Reading entry n − 1 of the list: O(1).

Over the whole doubling loop, by step 5 the sieve sizes sum to less than 2L, with L < 2pₙ (or L = 16 if the first sieve suffices). Total time is therefore O(pₙ log log pₙ), and space, the marker array of L + 1 entries plus the list of primes, is O(pₙ). By the prime number theorem pₙ is roughly n ln n, so in terms of n this is a little more than O(n log n). Compared with the brute force’s O(pₙ√pₙ), the factor √pₙ has shrunk to log log pₙ.

This post does not use the prime number theorem or any explicit bound on pₙ to justify the algorithm. The rough sizes above are only for estimating cost; correctness rests on the doubling alone.

Integer and index range.

• Stored values are at most the last L: marker-array indices run 0 to L, and every listed prime is at most L. Only transient loop values go past L: j after its final j += i, up to j + i < 2L, and the size L + 1 passed to the array. The doubling limit *= 2 runs only when a sieve falls short, so its largest result is the final L itself. For n = 10001 the final L is 131072 (see “Answer and checks”), so every value, transient ones included, stays below 2 × 131072 = 262144, far inside unsigned 32-bit range. The general condition is in the third bullet.
• The C++ code uses std::size_t for both values and indices. The marker array really occupies L + 1 entries, so on any platform where that array can be created, L is representable as a std::size_t.
• The product i × i is only computed after checking i ≤ L / i, so i × i ≤ L and it cannot overflow. In the crossing-out loop, j += i starts from j ≤ L and i ≤ L, so it reaches at most 2L; this, like L × 2 for the next sieve size, cannot overflow as long as L is at most half the maximum of std::size_t. An array that large runs out of memory first, so in practice memory is the real limit.
• nth_prime requires n ≥ 1, because it reads index n − 1 (see “Boundary cases” below).
• The standard specifies std::vector<bool> as a space-optimised specialisation, and common implementations store one entry per bit. Python’s bytearray uses one byte per entry. Python’s int has unlimited precision, so nothing overflows there.

C++17

#include <cstddef>
#include <iostream>
#include <vector>

// Sieve of Eratosthenes over 0..limit. Returns every prime <= limit in increasing order.
std::vector<std::size_t> primes_up_to(std::size_t limit) {
    std::vector<bool> composite(limit + 1, false);
    std::vector<std::size_t> primes;
    for (std::size_t i = 2; i <= limit; ++i) {
        if (composite[i]) {
            continue;
        }
        primes.push_back(i);
        // i <= limit / i is i * i <= limit without computing i * i.
        if (i <= limit / i) {
            for (std::size_t j = i * i; j <= limit; j += i) {
                composite[j] = true;
            }
        }
    }
    return primes;
}

// The n-th prime, for n >= 1. Doubles the sieve limit until the sieve holds at
// least n primes, so no bound on the n-th prime has to be known in advance.
std::size_t nth_prime(std::size_t n) {
    std::size_t limit = 16;
    while (true) {
        const std::vector<std::size_t> primes = primes_up_to(limit);
        if (primes.size() >= n) {
            return primes[n - 1];
        }
        limit *= 2;
    }
}

// Trial division. Used only to cross-check the sieve.
bool is_prime_by_trial(std::size_t m) {
    if (m < 2) {
        return false;
    }
    for (std::size_t d = 2; d <= m / d; ++d) {
        if (m % d == 0) {
            return false;
        }
    }
    return true;
}

int main() {
    if (nth_prime(1) != 2 || nth_prime(6) != 13) {
        return 1;
    }
    std::size_t count = 0;
    for (std::size_t m = 2; count < 100; ++m) {
        if (is_prime_by_trial(m)) {
            ++count;
            if (nth_prime(count) != m) {
                return 1;
            }
        }
    }
    std::cout << nth_prime(10001) << '\n';
    return 0;
}

is_prime_by_trial is the brute force’s trial division carried over unchanged, used only to check that the sieve gives the same answers for the first 100 primes. A mistake in the crossing-out range or in the doubling loop would be caught by this comparison. The small examples are checked with if, not assert. assert does nothing in a build where NDEBUG is defined, whereas the comparison and return 1 in an if are observable behaviour of the program, so -O2 optimisation cannot remove them.

Python 3

def primes_up_to(limit: int) -> list[int]:
    """Sieve of Eratosthenes over 0..limit. Returns every prime <= limit in increasing order."""
    composite = bytearray(limit + 1)
    primes = []
    for i in range(2, limit + 1):
        if composite[i]:
            continue
        primes.append(i)
        # The range is empty when i * i > limit.
        for j in range(i * i, limit + 1, i):
            composite[j] = 1
    return primes


def nth_prime(n: int) -> int:
    """The n-th prime, for n >= 1.

    Doubles the sieve limit until the sieve holds at least n primes, so no bound
    on the n-th prime has to be known in advance.
    """
    limit = 16
    while True:
        primes = primes_up_to(limit)
        if len(primes) >= n:
            return primes[n - 1]
        limit *= 2


def is_prime_by_trial(m: int) -> bool:
    """Trial division. Used only to cross-check the sieve."""
    if m < 2:
        return False
    d = 2
    while d <= m // d:
        if m % d == 0:
            return False
        d += 1
    return True


def main() -> None:
    if nth_prime(1) != 2 or nth_prime(6) != 13:
        raise SystemExit(1)
    count = 0
    m = 2
    while count < 100:
        if is_prime_by_trial(m):
            count += 1
            if nth_prime(count) != m:
                raise SystemExit(1)
        m += 1
    print(nth_prime(10001))


if __name__ == "__main__":
    main()

Same algorithm, same checks as the C++ version. There is no separate counterpart to the C++ if (i <= limit / i): when i × i > L, range(i * i, limit + 1, i) is empty, which has the same effect. bytearray(limit + 1) starts with every entry 0 (unmarked). If a small example fails, raise SystemExit(1) gives exit status 1. assert is not used because running with python -O removes it.

Answer and checks

The 10001st prime is 104743.

The parts you can follow by hand:

• Sieve size. The sieve grows 16, 32, 64, …, and since 2¹⁶ = 65536 < p₁₀₀₀₁ ≤ 2¹⁷ = 131072, sieves up to L = 65536 come up short and the loop stops at L = 131072. This agrees with L / 2 < pₙ ≤ L from step 5.
• Primality. The value above is odd, its digit sum is 19 so it is not a multiple of 3, and it does not end in 0 or 5 so it is not a multiple of 5. A complete check needs trial division by every number up to its square root (323 or less), and that is left to the programs rather than done by hand.

The checks come in three layers.

1. Small examples and the first 100 primes. Before printing the answer, both programs check p₁ = 2 and p₆ = 13, and check that the first 100 primes found by trial division equal the 1st through 100th primes produced by the sieve. If anything differs, they exit with a nonzero status. p₆ = 13 and p₇ = 17 were also traced by hand in the derivation.
2. Agreement between the two implementations. In this series’ repository a check script compiles the C++ code with g++ -std=c++17 -O2 -Wall -Wextra -Werror and runs it (warnings count as failures), and runs the Python code with python3 -I -B. Both outputs must match the stored answer.txt byte for byte, and a nonzero exit status or a timeout (1 second for C++, 10 seconds for Python) is a failure. By the series rules only code that passes this check is published.
3. The proof. That the sieve produces exactly the primes up to L, and that the doubling loop returns pₙ and terminates, was shown in the derivation and the correctness proof.

Boundary cases

• n = 1: returns 2, the first entry of the first sieve (L = 16). 1 is not prime and never enters the list, because the outer loop starts at i = 2. The small examples include p₁ = 2.
• n = 6: the first sieve holds exactly six primes, so the loop ends without doubling. That relies on the comparison being >=; writing > by mistake would sieve once more for nothing (with the same answer).
• n = 7: the first case that needs the sieve to grow; it returns 17 at L = 32.
• n = 0: there is no 0th prime, and the function requires n ≥ 1. Called with 0, the C++ n − 1 wraps around to the maximum unsigned value and the access goes out of range, while Python’s primes[-1] silently returns the last entry. This program only calls it with n ≥ 1, but reusing it elsewhere would need an input check first.
• Starting size: 16 was chosen only so that the problem’s example fits in the first sieve; any value of 1 or more is correct. Starting at 0 would never finish, since doubling 0 gives 0.
• Crossing out from i × i: starting at 2 × i gives the same answer but crosses out numbers that are already crossed out. Starting anywhere above i × i instead would leave composites behind and give wrong answers.

Going further

• Sieving once. With a proven explicit upper bound for pₙ you could sieve just once, but you would also have to check the range of n for which that bound actually holds. Doubling needs no such check and finishes safely at the price of at most a constant factor more work.
• Sieving odd numbers only. Handle 2 separately and store only odd numbers in the array, and memory and work drop to roughly half. Work out the mapping between index and number yourself (slot k holds 2k + 1).
• Not sieving again. Right now every doubling sieves from scratch. Could you instead use the primes already found to sieve only the new interval (L, 2L]? The key is where to start crossing out in each interval (a segmented sieve).
• Related problems. Problem 3 in this series found prime factors by trial division. Problem 10 deals with all primes below a fixed bound, so the sieve with a known bound applies directly.

Sources and license

• Project Euler, Problem 7 — the original problem. The “Problem” section of this post is my own summary of it and is shared under the same CC BY-NC-SA 4.0 license as the original. Attribution for the problem goes to Project Euler.
• Project Euler, Copyright — the license of the problem text and how to link to problems. Because of the non-commercial condition, this post carries no ads (ads: false).
• std::vector<bool> follows the C++ working draft [vector.bool], and assert and NDEBUG follow [assertions.assert]. For Python, assert and -O follow The assert statement, and bytearray follows Bytearray Objects.
• Mertens’ theorem and the prime number theorem are well-known results of number theory, quoted only to estimate the size of the cost; they are not proved here, and the correctness of the algorithm does not depend on either.
• The solution, proofs and code are my own. I did not consult Project Euler’s overviews, forum, member-only material or other solution sites. Motitory is not affiliated with or sponsored by Project Euler, and Project Euler has not reviewed this post.