More Prime-Counting Functions

July 26, 2011

Don’t be afraid of the Greek letters; you already know $\phi$ and $\pi$ , the $\sum$ is just a loop, the double $\sum \sum$ is just one loop nested within another, and the $\lfloor x \rfloor$ is just an indication that the arithmetic is done with integers rather than real numbers. Here’s our version of Meissel’s formula:

(define (meissel-pi n) (if (< n max-pi) (pi n) (let ((b (pi (isqrt n))) (c (pi (iroot 3 n)))) (let loop ((i (+ c 1)) (sum (+ (phi n c) (quotient (* (+ b c -2) (- b c -1)) 2)))) (if (< b i) sum (loop (+ i 1) (- sum (pi (quotient n (p i))))))))))

Lehmer’s formula isn’t much harder:

(define (lehmer-pi n) (if (< n max-pi) (pi n) (let ((a (pi (iroot 4 n))) (b (pi (isqrt n))) (c (pi (iroot 3 n)))) (let i-loop ((i (+ a 1)) (sum (+ (phi n a) (quotient (* (+ b a -2) (- b a -1)) 2)))) (if (< b i) sum (let* ((w (quotient n (p i))) (lim (pi (isqrt w))) (sum (- sum (pi w)))) (if (< c i) (i-loop (+ i 1) sum) (let j-loop ((j i) (sum sum)) (if (< lim j) (i-loop (+ i 1) sum) (j-loop (+ j 1) (- sum (pi (quotient w (p j))) (- j) 1)))))))))))

A library function to compute integer roots is given in the codepad source, which also includes all the code from the previous exercise. Our timing comparisons show that each successive function is about an order of magnitude faster than the previous one; in each case, we cleared the cache in the phi function, then computed each value in order:

legendre meissel lehmer (prime-pi n) 10^0 0.000 0.000 0.000 0 10^1 0.000 0.000 0.000 4 10^2 0.000 0.000 0.000 25 10^3 0.000 0.000 0.000 168 10^4 0.002 0.000 0.000 1229 10^5 0.010 0.002 0.002 9592 10^6 0.072 0.012 0.005 78498 10^7 0.739 0.059 0.021 664579 10^8 3.836 0.509 0.089 5761455 10^9 41.932 2.626 0.791 50847534 10^10 --- 12.579 2.163 455052511 10^11 --- 124.483 13.539 4118054813 10^12 --- --- 125.810 37607912018

The timings were done on a recent-vintage personal computer, running Petite Chez Scheme on Linux.

Like this:

Be the first to like this post.

Pages: 1 2

Posted by programmingpraxis

Filed in Exercises

3 Comments »

3 Responses to “More Prime-Counting Functions”

Graham said
July 26, 2011 at 1:26 PM
I think your write-up of Lehmer’s formula on the first page doesn’t match your code on the second. Specifically, what is b_i supposed to be? And is the fraction out front (b + a - 2) * (b - a + 1) / 2 or (b + c - 2) * (b - c + 1) / 2?
programmingpraxis said
July 26, 2011 at 1:53 PM
Oops. I didn’t do that very well, did I? It’s now fixed.

b_i = pi(isqrt(x/p_i) for a<i<=c

The fraction is (b+a-2)(b-a+1)/2
Graham said
July 26, 2011 at 2:38 PM
I’m completely off my game; Python is either giving me nonsense or complaining of exceeding the recursion depth. I’m going to start over, so I won’t have a solution for a bit.

Programming Praxis