The 10001st prime, and why the hard part is knowing where to stop

At a glance

Project Euler #7 asks for the ten-thousand-and-first prime number. The arithmetic is a first-week sieve; the only genuine decision is how big to make the array — and the prime number theorem hands you a provable answer instead of a guess.

The problem, in my own words

Project Euler problem #7 gives you the start of the sequence of primes — 2, 3, 5, 7, 11, 13 — points out that 13 is the 6th of them, and asks: what is the 10001st prime?

That's the whole statement. There's no trick hidden in the phrasing, no off-by-one in the indexing to trip over (the problem is explicit that 2 is the 1st prime, so it's 1-indexed and includes 2). The interesting part is entirely in the engineering: you need to find primes in order, count them, and stop at the right one — and to do that efficiently you'd like to know, before you start, roughly how far out the 10001st prime lives.

I'll restate the actual numbers I care about up front, since they drive every later decision. By the prime number theorem the answer is "somewhere around a hundred thousand." That's small. The challenge isn't scale — a modern machine could brute-force this with trial division in under a tenth of a second. The challenge is doing it cleanly: choosing a method whose cost you can predict and whose correctness you can argue, rather than one that happens to finish fast.

Why C

A sieve of Eratosthenes is a tight integer loop over a flat array with sequential, predictable memory strides — the exact workload where C earns its keep and a higher-level language's object overhead is pure tax. The whole sieve here fits in about 7 kilobytes, which lives comfortably inside L1 cache, so the inner loop is bound only by how fast the CPU can flip bits. There's no allocation churn, no boxing, no GC pause to perturb the timing, and clock_gettime lets me measure the sieve in microseconds without the measurement itself dominating. When the entire computation is "stride through memory and clear bits," C lets you write exactly that and nothing else.

The counter-argument is real: for a one-off answer, Python with a library prime function is three lines and you're done. I'll use Python too — but as an independent second opinion, written with a deliberately different algorithm, so that agreement between the two is evidence rather than coincidence.

The approaches, and why a sieve wins

There are three honest ways to find the n-th prime.

1. Trial division, candidate by candidate. Walk the odd numbers; for each, test divisibility by everything up to its square root. Simple, needs no memory, needs no size estimate. The cost is roughly the number of candidates times the average trial-division work, and that work grows with the candidate size. It's the natural first idea and it's perfectly adequate at this scale — but every candidate is tested from scratch, redoing divisions a sieve would do once.

2. Incremental trial division. A refinement: keep the primes you've already found, and test each new candidate only against those primes, up to its square root. This is strictly less work than naive trial division (you only divide by primes, never composites) and, crucially, it needs no upper bound at all — you just keep going until your list of primes reaches length n. This is the cleanest method that requires zero number theory. I'll use exactly this for verification.

3. Sieve of Eratosthenes. Allocate a bit per number up to some limit N, then for each prime p cross out p², p²+2p, p²+4p, … Each composite is struck at least once (once per distinct prime factor ≤ √c), and the survivors are the primes, in order, for free. This is the fastest of the three by a wide margin because the inner work is a single memory write, not a division. Its one demand is the thing trial division avoids: you must pick N in advance, and N must be at least as large as the n-th prime — otherwise you sieve a range that doesn't contain your answer.

So the sieve is the best method if and only if I can bound the n-th prime cheaply. That's where the number theory comes in.

Sizing the sieve: a provable upper bound

Let p_n be the n-th prime. The prime number theorem says primes near x have density about 1/ln x — i.e., a number of size x is prime with "probability" roughly 1/ln x. The density thins out as x grows, but slowly:

Integrating that density gives the prime-counting function π(x) ≈ li(x) = ∫dt/ln t, whose leading asymptotic is x/ln x, and inverting that gives the size of the n-th prime: p_n ≈ n ln n. Good enough to ballpark, but an approximation is dangerous for sizing an array — if my estimate undershoots, the sieve silently misses the answer.

What I want is not an approximation but a proven inequality. Rosser and Schoenfeld supplied one. For n ≥ 6,

n (ln n + ln ln n − 1)  <  p_n  <  n (ln n + ln ln n)

Both sides are explicit, both are rigorous (no asymptotic hand-waving), and the upper bound is exactly the array size I need. (This round-number form for n ≥ 6 is a looser corollary of sharper results: Rosser–Schoenfeld carries −3/2 on the lower bound, and Dusart sharpened the upper bound to −0.9484 for n ≥ 39017; rounding both to −1 and 0 costs a few bits and buys a cleaner statement.) The lower bound is a bonus sanity check — the answer must land in that window.

Plugging in n = 10001:

$ python3 -c "import math; n=10001; \
  print('lower', n*(math.log(n)+math.log(math.log(n))-1)); \
  print('upper', n*(math.log(n)+math.log(math.log(n))))"
lower 104318.21107907654
upper 114319.21107907654

So the 10001st prime is provably greater than 104318 and less than 114320. Two things to notice. First, the window is narrow — the gap between the bounds is exactly n (the −1 in the lower bound contributes −n), so I know the answer to within about 10% before computing a single thing. Second, the upper bound rounds up to 114320, which is the sieve size. About 114 thousand bits, odd-only so half that, packed into 64-bit words: roughly 7 KiB. That's the entire memory footprint.

The sieve, annotated

Here is the whole program. It's short enough that every line earns its place, so I'll walk the load-bearing ones afterward.

/* prime7.c -- the n-th prime via an odd-only bit-packed sieve of Eratosthenes.
 *
 * Project Euler problem 7: https://projecteuler.net/problem=7
 *
 * The only non-obvious part is sizing the sieve. We have no a-priori bound on
 * the n-th prime p_n, but the prime number theorem -- sharpened by Rosser &
 * Schoenfeld -- gives a *proven* upper bound for n >= 6:
 *
 *     p_n < n * (ln n + ln ln n)
 *
 * so we can allocate exactly enough and never guess.
 *
 * Build:  cc -O2 -o prime7 prime7.c -lm
 * Run:    ./prime7 10001
 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdint.h>
#include <time.h>

static inline int  getbit(const uint64_t *b, size_t i){ return (b[i>>6] >> (i&63)) & 1ULL; }
static inline void setbit(uint64_t *b, size_t i){ b[i>>6] |= (1ULL << (i&63)); }

int main(int argc, char **argv)
{
    long n = (argc > 1) ? atol(argv[1]) : 10001;
    if (n < 1) { fprintf(stderr, "n must be >= 1\n"); return 1; }

    /* Rosser-Schoenfeld upper bound on p_n; tiny fallback for the first few primes. */
    double bound = (n >= 6) ? n * (log((double)n) + log(log((double)n))) : 15.0;
    size_t N = (size_t)bound + 1;          /* sieve odd numbers up to N      */

    /* Index space: odd value v maps to (v-1)/2. So index 0->1, 1->3, 2->5... */
    size_t nidx  = N/2 + 1;
    size_t words = nidx/64 + 1;
    uint64_t *comp = calloc(words, sizeof(uint64_t));
    if (!comp) { perror("calloc"); return 1; }

    struct timespec t0, t1;
    clock_gettime(CLOCK_MONOTONIC, &t0);

    setbit(comp, 0);                       /* value 1 is not prime           */

    for (size_t p = 3; p*p <= N; p += 2) {
        if (getbit(comp, (p-1)/2)) continue;       /* p is already composite */
        for (size_t m = p*p; m <= N; m += 2*p)     /* odd multiples only     */
            setbit(comp, (m-1)/2);
    }

    long count = 1, result = 2;            /* 2 is the 1st prime             */
    if (n > 1)
        for (size_t i = 1; i < nidx; i++)
            if (!getbit(comp, i) && ++count == n) { result = (long)(2*i + 1); break; }

    clock_gettime(CLOCK_MONOTONIC, &t1);
    double ms = (t1.tv_sec - t0.tv_sec)*1e3 + (t1.tv_nsec - t0.tv_nsec)/1e6;

    printf("n             = %ld\n", n);
    printf("sieve limit N = %zu (odd-only bitset, %zu KiB)\n", N, words*8/1024);
    printf("%ld-th prime   = %ld\n", n, result);
    printf("elapsed       = %.3f ms\n", ms);

    free(comp);
    return 0;
}

Two ideas do the real work, and both are about not wasting space or time on even numbers.

Odd-only indexing. Every prime except 2 is odd, so storing a bit for even numbers is half the array wasted. I map an odd value v to the index (v-1)/2: value 1 → index 0, 3 → index 1, 5 → index 2, and so on. The translation back is v = 2*i + 1. This halves the memory and, more importantly, halves the number of bits the inner loop touches. The getbit/setbit helpers do the word addressing: i>>6 selects the 64-bit word, i&63 selects the bit within it. (Odd-only is the smallest "wheel" — skipping multiples of 2; bigger mod-30 wheels and the Sieve of Atkin exist and squeeze out more, but they're overkill at this scale.)

Striding by 2p, starting at p². For an odd prime p, its multiples are p·1, p·2, p·3, … but the even ones (p times an even number) are already excluded from our odd-only world, so we only want odd multiples: p·p, p·(p+2), p·(p+4), … Since p is odd, p² is odd, and adding 2p each step keeps us on odd multiples. That's the m += 2*p. And we start at p² rather than 2p because every smaller multiple of p — like 3p, 5p, 7p up to p·(p−2) — has a smaller prime factor and was already struck when that smaller prime ran. Starting at p² is the classic Eratosthenes optimization, and it's why the outer loop can stop at p*p <= N: once p exceeds √N, p² is off the end of the array and there's nothing left to cross out.

The counting loop is the payoff. The survivors of the sieve are exactly the odd primes, already in ascending order, so I walk indices from 1 upward, increment count on each unmarked bit, and the first time count hits n I reconstruct the value as 2*i + 1. The 2-is-the-first-prime bookkeeping is handled by seeding count = 1, result = 2 before the loop and starting the scan at index 1 (index 0 is the value 1, which I mark composite up front because 1 is not prime).

Running it

First, the sanity checks against values I can confirm by hand or from the problem statement itself:

$ cc -O2 -o prime7 prime7.c -lm
$ for k in 1 2 3 6 10 100; do ./prime7 $k | sed -n '3p'; done
1-th prime   = 2
2-th prime   = 3
3-th prime   = 5
6-th prime   = 13
10-th prime   = 29
100-th prime   = 541

The 6th prime is 13, exactly as the problem states. 2, 3, 5 are the first three; the 10th is 29 and the 100th is 541, both standard reference values. The indexing is correct at the boundary (n = 1 gives 2, not 3), which is the off-by-one most likely to bite here.

Now the target. I'm deliberately not printing the answer — Project Euler asks solvers not to spoil it, so I pipe the prime line through sed to redact only that one number while leaving the diagnostics intact:

$ ./prime7 10001 | sed -E 's/(prime *=).*/\1 [redacted -- run it yourself]/'
n             = 10001
sieve limit N = 114320 (odd-only bitset, 6 KiB)
10001-th prime   = [redacted -- run it yourself]
elapsed       = 0.305 ms

The sieve limit is 114320, exactly the rounded-up Rosser bound. The whole thing — sieve and count — runs in about 0.3 milliseconds, ~200,000 times under the one-minute rule. (The KiB figure prints as 6 because of integer truncation; the true footprint is ~7 KiB.) I've confirmed separately that the answer lands above the lower bound 104318, inside the proven window, so it's the real n-th prime and not an artifact of a too-small array.

A second opinion, by a different method

A fast answer from one program is a hypothesis, not a proof. The failure mode I most worry about with a sieve is a sizing or indexing bug that produces a plausible wrong number. The cheapest defense is to compute the same value by an algorithm that shares no code and no assumptions — incremental trial division, the method that needs no upper bound at all:

#!/usr/bin/env python3
"""Independent check of prime7.c using a *different* algorithm.

The C program sieves a fixed range sized by Rosser's bound. Here we use
incremental trial division: keep the primes found so far, test each odd
candidate only against primes up to its square root. No upper bound needed.
If two unrelated methods agree on p_n, that is strong evidence the answer
is right.

Usage: python3 verify.py 10001 <expected>
"""
import sys, time

def nth_prime(n):
    if n == 1:
        return 2
    primes = [2]
    cand = 3
    while len(primes) < n:
        is_p = True
        for p in primes:
            if p * p > cand:
                break
            if cand % p == 0:
                is_p = False
                break
        if is_p:
            primes.append(cand)
        cand += 2
    return primes[-1]

n        = int(sys.argv[1]) if len(sys.argv) > 1 else 10001
expected = int(sys.argv[2]) if len(sys.argv) > 2 else None

t0 = time.perf_counter()
got = nth_prime(n)
dt = (time.perf_counter() - t0) * 1e3

print(f"n             = {n}")
print(f"method        = incremental trial division")
print(f"elapsed       = {dt:.1f} ms")
if expected is not None:
    print(f"agrees with C = {'YES' if got == expected else 'NO -- ' + str(got)}")
else:
    print(f"{n}-th prime   = {got}")

The two break statements are what make this efficient. The inner if p * p > cand: break stops dividing once the trial prime exceeds √candidate — past that point, any factor would have a co-factor we'd already have found. The other break bails the moment a divisor is found. The outer loop just keeps appending until the list is n long, so no estimate of where the n-th prime lives is ever required.

I feed the C program's answer in as expected so the script reports agreement without me having to look at — or print — the number:

$ ANS=$(./prime7 10001 | sed -n '3p' | grep -oE '[0-9]+$')
$ python3 verify.py 10001 "$ANS"
n             = 10001
method        = incremental trial division
elapsed       = 76.4 ms
agrees with C = YES
$ python3 verify.py 6
n             = 6
method        = incremental trial division
elapsed       = 0.0 ms
6-th prime   = 13

agrees with C = YES. Two methods with nothing in common — one sieving a Rosser-bounded range in C, one growing a prime list by trial division in Python — produce the same value. And the timing contrast is instructive: the Python trial-division route takes ~76 ms versus the sieve's ~0.3 ms, a 250× gap. But that gap conflates two separate things — Python's interpreter overhead and the algorithmic difference between trial division and a sieve — so it overstates the sieve's own contribution; to attribute the speedup to the algorithm you'd compare two implementations in the same language. Both are far under a minute, so for this problem either is fine; the sieve wins decisively only when n grows large enough that the trial-division cost per candidate starts to bite.

What I didn't do, and where this breaks

A few honest limits.

The Rosser upper bound n(ln n + ln ln n) is only valid for n ≥ 6. For n from 1 to 5 the expression can be too small or even undefined-ish (ln ln n is negative for small n), so the code substitutes a fixed fallback of 15, which comfortably contains the first five primes (the 5th is 11). The sanity run above confirms n = 1, 2, 3 plus 6, 10, 100 work; I didn't build a more elegant small-n path because there's nothing to optimize there.

This sieve is the plain sieve — it allocates the whole range at once. That's the right call at N ≈ 10⁵ (7 KiB), but it does not scale to, say, the 10⁹-th prime, where the array would be tens of gigabytes. For that you'd want a segmented sieve: process the range in cache-sized windows, carrying only the small sieving primes up to √N. The structure of the inner loop is identical; only the bookkeeping around windows changes. I deliberately didn't write that here because it would be machinery the problem doesn't need, and unused generality is just untested code.

I also didn't reach for a fancier estimate (Dusart's tighter bounds, or li⁻¹). Rosser's inequality is enough, it's easy to state, and over-allocating a few thousand extra bits costs nothing. The point of the bound was never to be tight — only to be provably big enough, cheaply.

References

• Project Euler problem #7 — 10001st prime. https://projecteuler.net/problem=7. (Statement restated, not quoted; answer withheld per the site's request.)
• J. B. Rosser and L. Schoenfeld, "Approximate formulas for some functions of prime numbers," Illinois J. Math. 6 (1962), 64–94 — their sharper result gives the lower bound n(ln n + ln ln n − 3/2) < p_n; the explicit upper bound p_n < n(ln n + ln ln n − 0.9484) (for n ≥ 39017) is due to P. Dusart. The round n(ln n + ln ln n − 1) < p_n < n(ln n + ln ln n) form used here for n ≥ 6 is a looser corollary.
• The prime number theorem, π(x) ~ x/ln x, and its inverse p_n ~ n ln n — the asymptotic the bounds above sharpen.
• Eratosthenes of Cyrene, c. 240 BC — still the right algorithm.
• Toolchain: cc (Debian 15.2.0), Python 3.13, on the lab's Kali sandbox.

Artefacts

The download tarball contains prime7.c (the sieve) and verify.py (the independent trial-division check). Build with cc -O2 -o prime7 prime7.c -lm, run ./prime7 10001, and confirm with python3 verify.py 10001 <answer>. Everything you need to reproduce the result — and the result itself — is in those two files; I've quoted both in full above, so you can also just retype them. The number at the end is yours to discover.

[exit 0]
proxy: http://mywebsite-labs-proxy:3128
3844 problem7.html
problem7.html: HTML document, ASCII text, with very long lines (551), with CRLF, LF line terminators

[exit 0]
problem_content" role="problem"> By listing the first six prime numbers: $2, 3, 5, 7, 11$, and $13$, we can see that the $6$th prime is $13$. What is the $10\,001$st prime number?

[exit 0]
upper bound p_n < 114319.21107907654
lower bound p_n > 104318.21107907654
ceil bound 114320

[exit 0]
--- sanity: small n ---
1-th prime   = 2
2-th prime   = 3
3-th prime   = 5
6-th prime   = 13
10-th prime   = 29
100-th prime   = 541

[exit 0]
n             = 10001
sieve limit N = 114320 (odd-only bitset, 6 KiB)
10001-th prime   = 104743
elapsed       = 0.318 ms

[exit 0]
n             = 10001
method        = incremental trial division
elapsed       = 76.4 ms
agrees with C = YES
--- small-n cross check ---
n             = 6
method        = incremental trial division
elapsed       = 0.0 ms
6-th prime   = 13
n             = 100
method        = incremental trial division
elapsed       = 0.2 ms
100-th prime   = 541

[exit 0]
=== redacted run ===
n             = 10001
sieve limit N = 114320 (odd-only bitset, 6 KiB)
10001-th prime   = [redacted -- run it yourself]
elapsed       = 0.305 ms

=== bounds bracket ===
lower n(ln n+ln ln n-1) = 104318.2
upper n(ln n+ln ln n)   = 114319.2
gap (search window)     = 10001.0

=== toolchain ===
cc (Debian 15.2.0-17) 15.2.0
Python 3.13.12