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 `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:

(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.

(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.