Saturday, February 16, 2008

Doing Eratosthenes Justice

Yesterday I showed how to use some closures and macros to implement (very inefficient) lazy lists in Lisp. One of my examples produced an infinite list of prime numbers.

BlueStorm quite correctly pointed out that it was O(n^2), while the classic sieve of Eratosthenes is O(n log n). Always up for a challenge, I went ahead wrote a proper sieve.
(def merge_ (xs ys)
(with (x (car_ xs) y (car_ ys))
(if (< x y) (cons_ x (merge_ (cdr_ xs) ys))
(> x y) (cons_ y (merge_ xs (cdr_ ys)))
(cons_ x (merge_ (cdr_ xs) (cdr_ ys))))))

(def diff_ (xs ys)
(with (x (car_ xs) y (car_ ys))
(if (< x y) (cons_ x (diff_ (cdr_ xs) ys))
(> x y) (diff_ xs (cdr_ ys))
(diff_ (cdr_ xs) (cdr_ ys)))))

(def merge-multiples (li)
(with (xs (car_ li) rest (cdr_ li))
(cons_ (car_ xs)
(merge_ (cdr_ xs) (merge-multiples rest)))))

(def multiples (p) (map_ [* _ p] (count-from p)))

(= primes
(cons_ 2 (cons_ 3 (cons_ 5 (diff_ (count-from 7) non-primes)))))

(= non-primes (merge-multiples (map_ multiples primes)))
A typical optimization is to ignore multiples of 2, as follows.
(def count-from-by (n by) (cons_ n (count-from-by (+ n by) by)))

(def multiples (p) (map_ [* _ p] (count-from-by p 2)))

(= primes
(cons_ 2 (cons_ 3 (cons_ 5 (diff_ (count-from-by 7 2) non-primes)))))

(= non-primes (merge-multiples (map_ multiples (cdr_ primes))))

1 comment:

Michael Finney said...

Cool!

Wikipedia - Sieve of Eratosthenes has a colorful explanation. It's worth a peek too.