The Nth Prime
August 2, 2011
To calculate the nth prime we bracket n in loop1, use binary search in loop2, then calculate the prime-counting function and count up or down to the requested n in loop3:
(define (nth-prime n)
; (define (approx-pi x) (inexact->exact (round (riemann-r x))))
(define (approx-pi x) (inexact->exact (round
(- (logint x) (/ (logint (sqrt x)) 2) (/ (logint (expt x 1/3)) 3)))))
(if (< n 5) (vector-ref (vector 2 3 5 7) (- n 1))
(let loop1 ((lo 1) (hi 2))
(if (< (approx-pi hi) n) (loop1 hi (* hi 2))
(let loop2 ((lo lo) (hi hi))
(let* ((mid (quotient (+ lo hi) 2)) (mid-pi (approx-pi mid)))
(cond ((< n mid-pi) (loop2 lo mid))
((< mid-pi n) (loop2 mid hi))
(else (let* ((m (if (prime? mid) mid (prev-prime mid))) (k (prime-pi m)))
(cond ((< k n) (let loop3 ((k k) (m m))
(if (= n k) m (loop3 (+ k 1) (next-prime m)))))
((< n k) (let loop3 ((k k) (m m))
(if (= k n) m (loop3 (- k 1) (prev-prime m)))))
(else m)))))))))))
The test for n less than 5 is needed because the binary search “steps over” the correct values when n is so small. The call to (prev-prime mid) is necessary because the result of approx-pi is not necessarily prime.
This function makes a very good example of the trade-offs necessary in developing software. The “right” method uses the Riemann function to approximate π, as shown in the commented-out line of code. But that function is called many times in loop1 and loop2, and it’s too slow. A better solution uses the weaker approximation of the first three terms of the Riemann function at the cost of more calls to the next-prime and prev-prime functions but still provides a much faster function.
The millionth prime is 15485863; the billionth prime is 22801763489:
> (nth-prime #e1e6)
15485863
> (nth-prime #e1e9)
22801763489
We used all the machinery of the three previous exercises plus the Baillie-Wagstaff primality checker. You can run the program at http://programmingpraxis.codepad.org/KgVfGvEG.
Programming Praxis 
[…] primes, not knowing the largest prime in the set. One method is to calculate the nth prime, using a prime-counting formula, then use the Sieve of Eratosthenes to generate the list. Or instead of calculating the nth prime […]
I just came across this, and I think it’s a great idea and can be done pretty easily with just simple functions (a good primality tester being one of them, which as we’ve seen through many of your exercises is a really handy tool tool to have). I will point out that the GNU MPFR library (available in C, Perl, and Python to name a few) has a fast high-precision Riemann Zeta function, leaving only the relatively straightforward Riemann R calculation to do.
Rather than the search on Riemann R, an easier but somewhat slower way would be to use Cipolla’s nth prime approximation (Dusart 1999 page 12, among other sources) using, for instance, 3 terms (m=2). Then do Lehmer’s method on that and count up/down. Simple and straightforward, but we don’t get nearly as close as as using Riemann’s R function. I suppose it depends on how simple you want the code, how long Riemann R takes to compute, and how fast we can count.
My original solution, before doing the Lehmer prime counting function, was just optimizing the straightforward approach — fast segment sieve and fast bit counting. Once I added Lehmer’s method, given I had a fast segmented prime count, the method I used: compute a lower limit (Dusart 2010, page 2), do Lehmer’s method to get that count, then do the sieve from that point. The Dusart lower limit is cheap and only has to be computed once. Of course it doesn’t get as close as a search on Riemann’s R. The other advantage is that segment counting is pretty fast vs. prev_prime / next_prime.
Looking at 10^12 as an example:
– Using Riemann R, we come up with 29996225393466. Lehmer’s count gives us 1000000035884, so we call prev_prime 35885 times and get 29996224275833. All but 0.5s is spent doing the prime count (I just used a binary search). Total time in my code 11.8s.
– Using Cipolla’s approximation, we immediately get 29996300456744. Lehmer’s count gives us 1000002454213, then we count down to 29996224275833. Total time: 43.6s.
– Using the lower+sieve method, we immediately get 29994075087860. Lehmer’s count gives us 999930740354. Since we used the lower limit, we know this will always be <= the target. We sieve from there and get 29996224275833. Total time: 19.2s.
– Using just the sieve, we churn away until we get 29996224275833. Primesieve on 12 cores takes over 33 minutes, which goes to show just how much faster Lehmer's method is (especially as threaded primesieve is super fast as far as sieves go).
It looks like, with my implementation, things are a wash in performance between the R-search+prev_prime/next_prime and lower+sieve methods until about 10^11, where your Riemann R method starts coming out ahead. For 10^13th it's over 3x faster. Definitely something I'll want to use at some point.