## Counting Primes

### February 9, 2016

Here is our translation of the Hardy-Wright function into Scheme, accumulating the factorials as it goes:

```(define (prime-pi n)
(if (< n 2) 0 (if (< n 3) 1 (if (< n 4) 2
(let loop ((j 3) (j2fact 1) (pi -1))
(if (< n j) pi
(let ((j2fact (* j2fact (- j 2))))
(loop (+ j 1) j2fact (+ pi (- j2fact
(* j (quotient j2fact j))))))))))))```

That works properly, but it’s extremely slow even when n is only moderately large because the factorials quickly become huge; for instance, `(prime-pi 1000000)`, which isn’t really a very big number, takes over three hours on my computer at home:

```> (time (prime-pi 1000000))    24330 collections
11779966 ms elapsed cpu time, including 20871 ms collecting
11788363 ms elapsed real time, including 20896 ms collecting
3331524960608 bytes allocated, including 3331252634800 bytes reclaimed
78498```

By comparison, enumerating those primes with a sieve takes less than a second.

You can run the program at http://ideone.com/fJJtwD.

Pages: 1 2

### 7 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)))))
```
7. […] this weekend I was going through the Programming Praxis website, and this problem caught my […]