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() memory when O() 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 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 were first. I made another tweak to accumulate the first primes into a product so that I could check all of those against with a
(gcd n product-of-first-k) where was constrained so that the product stayed below
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
Ignoring the special-case at the beginning for two, my current algorithm looks like this:
(let ((nixed (make-instance 'cl-heap:fibonacci-heap
(loop :do (incf guess 2) :until (%prime-p guess nixed))
(cl-heap:add-to-heap nixed (cons (* guess guess) (* guess 2)))
%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.
(let ((prime t)
(loop :do (setf n (cl-heap:pop-heap nixed))
:while (and n (= (car n) guess))
(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)))
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.