Everybody loves the Sieve of Eratosthenes

unsorted — cgrand, 30 July 2009 @ 11 h 27 min

If I judge by traffic logs for this blog, many newcomers want to compute prime numbers in Clojure.

A recent thread on the mailing list prompted me to test an idea I had for a while: to use a map to implement the Sieve of Eratosthenes.

The first implementation was this one:

(defn primes [max]
  (let [enqueue (fn [sieve n factor]
                  (let [m (+ n factor)]
                    (assoc sieve m
                      (conj (sieve m) factor))))
        next-sieve (fn [sieve candidate]
                     (if-let [factors (sieve candidate)]
                       (reduce #(enqueue %1 candidate %2)
                         (dissoc sieve candidate)
                         factors)
                       (enqueue sieve candidate candidate)))]
    (apply concat (vals (reduce next-sieve {} (range 2 max))))))

where the sieve is a map from the next non-prime numbers to their factors. It’s so naive that even numbers are tested for primality but it doesn’t perform that bad: on my box, it takes 3s to compute all primes below 1,000,000 (it’s roughly as fast as clojure.contrib.lazy-seqs/primes when the seq isn’t yet memoized).

I wasn’t happy with the way I handled the case where a non-prime was already checked off (ie was a key of the map): I was conjing onto the list of prime factors for this number. If instead I tried to check off n+p (where n is the already known non-prime and p the current prime) or n+2p or n+3p… until I found a yet unknown non-prime I wouldn’t need to maintain a list of factors, nor to conj or reduce over it. And I was hoping that less allocations would yield better perfs.

Here is the second iteration:

(defn primes2 [max]
  (let [enqueue (fn [sieve n factor]
                  (let [m (+ n factor)]
                    (if (sieve m)
                      (recur sieve m factor)
                      (assoc sieve m factor))))
        next-sieve (fn [sieve candidate]
                     (if-let [factor (sieve candidate)]
                       (-> sieve
                         (dissoc candidate)
                         (enqueue candidate factor))
                       (enqueue sieve candidate candidate)))]
    (vals (reduce next-sieve {} (range 2 max)))))

and it computes all the primes below 1,000,000 in 1.8s instead of 3s but it still tests even numbers for primality.

primes3 is primes2 modified to only test odd numbers:

(defn primes3 [max]
  (let [enqueue (fn [sieve n factor]
                  (let [m (+ n (+ factor factor))]
                    (if (sieve m)
                      (recur sieve m factor)
                      (assoc sieve m factor))))
        next-sieve (fn [sieve candidate]
                     (if-let [factor (sieve candidate)]
                       (-> sieve
                         (dissoc candidate)
                         (enqueue candidate factor))
                       (enqueue sieve candidate candidate)))]
    (cons 2 (vals (reduce next-sieve {} (range 3 max 2))))))

and it computes the same list of primes in 1.5s.

Out of curiosity, I wrote a lazy version of primes3:

(defn lazy-primes3 []
  (letfn [(enqueue [sieve n step]
            (let [m (+ n step)]
              (if (sieve m)
                (recur sieve m step)
                (assoc sieve m step))))
          (next-sieve [sieve candidate]
            (if-let [step (sieve candidate)]
              (-> sieve
                (dissoc candidate)
                (enqueue candidate step))
              (enqueue sieve candidate (+ candidate candidate))))
          (next-primes [sieve candidate]
            (if (sieve candidate)
              (recur (next-sieve sieve candidate) (+ candidate 2))
              (cons candidate 
                (lazy-seq (next-primes (next-sieve sieve candidate) 
                            (+ candidate 2))))))]
    (cons 2 (lazy-seq (next-primes {} 3)))))

and, surprisingly (better locality?), it computes the primes below 1,000,000 in 1s.

9 Comments »

  1. Here's my naive, ugly and untested take on the problem:
    (defn not-divisible-by?[num denum]
    (not (= (mod num denum) 0)))

    (defn primes
    ([pos] (primes (vec (range 2 pos)) 0))
    ([clist pos]
    (cond
    (>= pos (count clist)) clist
    :else
    (let [val (nth clist pos)
    starting (inc pos)]
    (recur (into (subvec clist 0 starting) (vec (filter #(not-divisible-by? % val) (subvec clist starting (count clist))))) starting)))))

    Comment by zemariamm — 15 September 2009 @ 4 h 36 min
  2. [...] reading this post [...]

  3. [...] favorite solution was by Cristophe Grande, not just because he followed a similar approach to me, but also because he was very clever in his [...]

  4. [...] of SoE. I would recommend checking out Christophe Grand’s treatise on the subject titled Everybody loves the Sieve of Eratosthenes for a great discussion on writing real world prime sieves in [...]

  5. [...] Clojure is very big on laziness… and I’m feeling to lazy to re-implement the Sieve in a one-liner. I’m sure it could be done! In the meantime, I’ll point you to cgrand’s great article about the Sieve of Eratosthenes in clojure. [...]

  6. thanks for this good article..

    Comment by heiji — 27 August 2011 @ 16 h 15 min
  7. [...] I wasn’t happy with this performance, especially after seeing Christophe Grand’s solution and timing on his blog post Everybody loves the Sieve of Eratosthenes. [...]

  8. アバクロ 通販

    Comment by アバクロ パンツ — 11 September 2013 @ 12 h 06 min
  9. Web ? Le compte aidé m’une affaire acceptable. J’avais été petit le peu
    a mis au courant de ceci votre émission offert brillant clair la concept d’en haut|salut de Hello d’en haut|salut de MON et l’époux de
    JE aime absolument votre blog et trouve beaucoup de beaucoup de
    la plupart du de votre poste pour être précisément ce que je cherche.

    peut vous Les écrivains d’invité d’offre pour écrire le contenu pour vous ?
    Je n’aurais pas des objections publishing une poste ou élaborant sur
    certains de du expose vous écrit dans les égards d’à ici.
    Encore, impressionnant website! de Howdy !
    Utilisez-vous le Gazouillement ? J’aimerais vous suivre si cela serait approuve.

    Je suis definitely appréciant votre blog et attend avec impatience nouveau
    updates.

    Comment by hay day hack — 5 October 2014 @ 2 h 01 min

RSS feed for comments on this post. TrackBack URI

Leave a comment

(c) 2014 Clojure and me | powered by WordPress with Barecity