Faster Primes with Heaps January 28th, 2013
Patrick Stein

I read a nice blog post this morning about implementing the Sieve of Eratosthenes in Lisp. It got me thinking about something that I did a few weeks back.

A few weeks back, I needed to generate all of the primes up to a given number. The number was big enough that I didn’t want to keep an array of bytes (or even bits) big enough to do the sieving. I just had an aversion to keeping around O(n) memory when O(n / \ln n) would suffice. I also had a hankering (probably because of Let Over Lambda) to implement it as a generator—a function that each time you called it gave you the next prime.

My first version amounted to trial-division against the entire list of primes generated so far. For the number of primes that I wanted to generate, this took more than nine minutes. Unacceptable.

I made some tweaks to it to ensure that I was only checking n against the previously generated primes that were less than or equal to (isqrt n). I made some other tweaks to keep my previously generated primes sorted so the most-likely to be a divisor of n were first. I made another tweak to accumulate the first k primes into a product so that I could check all of those against n with a (gcd n product-of-first-k) where k was constrained so that the product stayed below most-positive-fixnum.

Those tweaks got me down to around two minutes. That was still more than twice as long as I wanted it to take.

Now, with the help of CL-HEAP and the Rosetta-Code-approved optimization mentioned in that blog post, I have it down to thirty seconds with no trial division, no square roots, and only one multiplication per prime.

I keep a heap of cons cells. The car is the next number that will be crossed out by this prime and the cdr is how much to increment the car once I’ve crossed out the car.

Ignoring the special-case at the beginning for two, my current algorithm looks like this:

(defun prime-generator ()
  (let ((nixed (make-instance 'cl-heap:fibonacci-heap
                              :key #'first
                              :sort-fun #'<))
        (guess 1))
   
    (lambda ()
      (loop :do (incf guess 2) :until (%prime-p guess nixed))
      (cl-heap:add-to-heap nixed (cons (* guess guess) (* guess 2)))
      guess)))

The %prime-p function pulls cons cells off of the heap until it finds whose car isn’t equal to the current guess. For each of them where the car was equal, it increments the car by its cdr. Then, it puts all of the cons cells it removed back on the heap. The guess was prime if it wasn’t equal to any of the numbers on the heap.

(defun %prime-p (guess nixed)
  (let ((prime t)
        (n nil))
    (loop :do (setf n (cl-heap:pop-heap nixed))
       :while (and n (= (car n) guess))
       :do (progn
             (incf (car n) (cdr n))
             (cl-heap:add-to-heap nixed n)
             (setf prime nil))
       :finally (when n (cl-heap:add-to-heap nixed n)))
    prime))

I could probably save some time and cons-ing by peeking at the top of the heap rather than always gathering or countless other ways, but this more than exceeded my expectations for today.

Updates In Email

Email:

RSS Feeds

l