Counting Primes

February 9, 2016

We have studied functions for counting the prime numbers less than a given number in two previous exercises. All of them were based on Legendre’s phi function that counts the numbers less than x not stricken by sieving with the first a primes.

Today we look at a rather different method of counting the primes that is due to G. H. Hardy and Edward M. Wright. Their method is based on factorials, and their formula is

\pi(n) = -1 + \sum_{j=3}^{n} \left( (j-2)! - j\lfloor \frac{(j-2)!}{j} \rfloor \right)

for n > 3, where ⌊n⌋ is the greatest integer less than n. The expression inside the big parentheses is 1 when n is prime and 0 when n is composite, by Wilson’s Theorem, which states that an integer n > 1 is prime if and only if (n − 1)! ≡ −1 (mod n); the theorem was first stated by Ibn al-Haytham (c. 1000AD), but is named for John Wilson, who first published it in 1770, and it was first proved by Lagrange in 1771.

Your task is to write a program that counts primes by the Hardy-Wright method. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

6 Responses to “Counting Primes”

  1. Very interesting algo, albeit somewhat impractical. It was a good refresher though for the BigInt type (Julia language) that makes it usable for inputs over 23.

  2. Jan Van lent said
    (use-package :iterate)
    (defun prime-count (n)
      (+ -1 (iter (for j from 3 to n)
                  (for f first 1 then (* f (- j 2)))
                  (summing (- f (* j (floor f j)))))))
  3. Jan Van lent said

    My previous comment has a solution in Common Lisp.

  4. Mike said
    from itertools import accumulate, count, islice
    from operator import mul
    def pi(n):
        factorials = accumulate(count(1), mul)
        summands = (f - (f//j)*j for j,f in enumerate(factorials, 3))
        return sum(islice(summands, n-2)) - 1
  5. DanN said
    import functools
    def fact(n):
        return functools.reduce(lambda x, y: x * y, list(range(1, n + 1)))
    def primeCounting(n):
        return (-1 + sum([(fact(j - 2) - j * (fact(j - 2) // j)) for j in range(3, n + 1)]))
  6. Jussi Piitulainen said

    Two Scheme implementations, or one in two ways. I felt quite stupid for quite a while before I understood that I want to update the factorial with j-1 instead of j-2 when I also update j. (Now I don’t understand why the model implementation also works :)

    (define (pi n) ;for n > 3
      (do ((S 0 (+ S (- f (* j (quotient f j)))))
           (f 1 (* (- j 1) f)) ; sic!
           (j 3 (+ j 1)))
          ((< n j) (- S 1))))
    (define (pi1 n) ;for n > 3
      (let loop ((S 0) (f 1) (j 3))
        (if (< n j)
            (- S 1)
            (loop (+ S (- f (* j (quotient f j))))
                  (* (- j 1) f) ; sic!
                  (+ j 1)))))

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: