A story about a simple problem that wasn't simple, and a bug that cost me 2 minutes of runtime.
The problem
Given N (1 ≤ N ≤ 10^9), print the N-th prime number.
Samples: N = 3 gives 5, N = 30 gives 113, N = 300000000 gives 6461335109.
Looks like a textbook sieve problem, right? Let's see.
Attempt 1: the sieve
The N-th prime is about N · ln(N). For N = 10^9 the answer is 22801763489.
A sieve has to cover all numbers up to that, roughly 2.3 · 10^10. A plain bool array would need 23 GB of memory, so it's dead on arrival.
Attempt 2: a better sieve
The usual tricks help a lot:
- store odd numbers only (half the memory and work);
- use a segmented sieve with ~32 KB blocks that fit in the L1 cache;
- start crossing out from p*p.
Memory is no longer a problem, but the work is: about 10^10 elementary operations is still tens of seconds. Wheel factorization and bitsets give constant-factor gains, but we need something like a 500x speedup. Time to change the algorithm.
The real idea: don't sieve, count!
Instead of generating all primes, we ask a different question: where is the N-th prime, and how many primes are below that point? The solution has three steps.
Step 1. Guess where the answer is
Let π(x) be the number of primes ≤ x, and let li(x) be the logarithmic integral:
li(x) = integral from 0 to x of dt / ln(t)
We know π(x) ≈ li(x) with a very small error (around √x · ln(x) in practice, far less than x). So the N-th prime is close to the solution of li(x) = N. Solve it with Newton's method, using li'(x) = 1 / ln(x):
x = x - (li(x) - N) * ln(x)
To evaluate li(x) I used Ramanujan's series:
li(x) = γ + ln(ln(x)) + √x * Σ (n = 1 to ∞) of a(n) * b(n)
a(n) = (-1)^(n-1) * (ln x)^n / (n! * 2^(n-1))
b(n) = 1/1 + 1/3 + 1/5 + ... + 1/(2k+1), where k goes up to floor((n-1)/2)
It's an alternating series with large intermediate terms, so use long double.
Step 2. Count primes exactly with the Lucy_Hedgehog method
Now we need the exact value of π(x) for x ≈ 2 · 10^10. There is a beautiful O(x^(3/4)) algorithm for it.
Let S(v, p) be the number of integers in [2, v] that are either prime or have all prime factors greater than p. Initially S(v, 1) = v — 1. When we process a prime p:
S(v, p) = S(v, p-1) - ( S(floor(v/p), p-1) - S(p-1, p-1) )
for all v ≥ p*p. The key observation is that we only ever need v of the form floor(x/i), and there are only about 2√x such values. So the arrays have size about 1.5 · 10^5 and everything fits in cache. After processing all p ≤ √x, we get S(x, √x) = π(x).
Step 3. Fix the small error
Now we know c = π(X) for our guess X. If c < N, sieve a small window to the right of X and walk until we meet the N-th prime. If c ≥ N, sieve to the left and walk back c — N primes. The window is only a couple of million numbers, which is trivial.
Results
- N = 1: answer 2, time ~0.01 s
- N = 30: answer 113, time ~0.03 s
- N = 10^6: answer 15485863, time ~0.09 s
- N = 3 · 10^8: answer 6461335109, time ~0.06 s
- N = 10^9: answer 22801763489, time ~0.11 s
Compared to tens of seconds for the best sieve, that's a speedup of about three orders of magnitude.
The bug that cost me 2 minutes
My first version was correct but terribly slow: 33 seconds for N = 3 · 10^8 and 119 seconds for N = 10^9. The answers were right, so the algorithm was fine and something else was wrong.
The culprit: in the li(x) series I initialized the coefficient as L/2 instead of L (the n = 1 term is L^1 / (1! * 2^0) = L, where L = ln x), so li(x) was off by a factor of two. Newton's method converged to the wrong place (about 1.3 · 10^10 instead of 6.5 · 10^9 for N = 3 · 10^8). The final window walk then had to cross billions of numbers, and the "small fix" step turned into a full sieve.
Lesson: when the answer is correct but the time is not, check that your estimate is actually close to the truth. A cheap sanity check, comparing the guess with the known value, would have found this in a minute.
Code
Takeaways
- If a sieve needs more than about 10^9 numbers, think about counting instead of generating.
- li(x) is an excellent approximation of π(x), so its inverse gives a precise guess for the N-th prime.
- Lucy_Hedgehog computes π(x) in O(x^(3/4)) with tiny memory and about 20 lines of code.
- Always sanity-check intermediate estimates, not just the final answers.
Thanks for reading! Questions and improvements are welcome in the comments. For example, can anyone push this further with Meissel-Lehmer or a faster π(x) in O(x^(2/3))?
