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

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:

WordPress.com Logo

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

Google photo

You are commenting using your Google 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 )

Connecting to %s

%d bloggers like this: