Fast Prime Number Generation in Clojure

Viewed 17003

I've been working on solving Project Euler problems in Clojure to get better, and I've already run into prime number generation a couple of times. My problem is that it is just taking way too long. I was hoping someone could help me find an efficient way to do this in a Clojure-y way.

When I fist did this, I brute-forced it. That was easy to do. But calculating 10001 prime numbers took 2 minutes this way on a Xeon 2.33GHz, too long for the rules, and too long in general. Here was the algorithm:

(defn next-prime-slow
    "Find the next prime number, checking against our already existing list"
    ([sofar guess]
        (if (not-any? #(zero? (mod guess %)) sofar)
            guess                         ; Then we have a prime
            (recur sofar (+ guess 2)))))  ; Try again                               

(defn find-primes-slow
    "Finds prime numbers, slowly"
    ([]
        (find-primes-slow 10001 [2 3]))   ; How many we need, initial prime seeds
    ([needed sofar]
        (if (<= needed (count sofar))
            sofar                         ; Found enough, we're done
            (recur needed (concat sofar [(next-prime-slow sofar (last sofar))])))))

By replacing next-prime-slow with a newer routine that took some additional rules into account (like the 6n +/- 1 property) I was able to speed things up to about 70 seconds.

Next I tried making a sieve of Eratosthenes in pure Clojure. I don't think I got all the bugs out, but I gave up because it was simply way too slow (even worse than the above, I think).

(defn clean-sieve
    "Clean the sieve of what we know isn't prime based"
    [seeds-left sieve]
    (if (zero? (count seeds-left))
        sieve              ; Nothing left to filter the list against
        (recur
            (rest seeds-left)    ; The numbers we haven't checked against
            (filter #(> (mod % (first seeds-left)) 0) sieve)))) ; Filter out multiples

(defn self-clean-sieve  ; This seems to be REALLY slow
    "Remove the stuff in the sieve that isn't prime based on it's self"
    ([sieve]
        (self-clean-sieve (rest sieve) (take 1 sieve)))
    ([sieve clean]
        (if (zero? (count sieve))
            clean
            (let [cleaned (filter #(> (mod % (last clean)) 0) sieve)]
                (recur (rest cleaned) (into clean [(first cleaned)]))))))

(defn find-primes
    "Finds prime numbers, hopefully faster"
    ([]
        (find-primes 10001 [2]))
    ([needed seeds]
        (if (>= (count seeds) needed)
            seeds        ; We have enough
            (recur       ; Recalculate
                needed
                (into
                    seeds    ; Stuff we've already found
                    (let [start (last seeds)
                            end-range (+ start 150000)]   ; NOTE HERE
                        (reverse                                                
                            (self-clean-sieve
                            (clean-sieve seeds (range (inc start) end-range))))))))))

This is bad. It also causes stack overflows if the number 150000 is smaller. This despite the fact I'm using recur. That may be my fault.

Next I tried a sieve, using Java methods on a Java ArrayList. That took quite a bit of time, and memory.

My latest attempt is a sieve using a Clojure hash-map, inserting all the numbers in the sieve then dissoc'ing numbers that aren't prime. At the end, it takes the key list, which are the prime numbers it found. It takes about 10-12 seconds to find 10000 prime numbers. I'm not sure it's fully debugged yet. It's recursive too (using recur and loop), since I'm trying to be Lispy.

So with these kind of problems, problem 10 (sum up all primes under 2000000) is killing me. My fastest code came up with the right answer, but it took 105 seconds to do it, and needed quite a bit of memory (I gave it 512 MB just so I wouldn't have to fuss with it). My other algorithms take so long I always ended up stopping them first.

I could use a sieve to calculate that many primes in Java or C quite fast and without using so much memory. I know I must be missing something in my Clojure/Lisp style that's causing the problem.

Is there something I'm doing really wrong? Is Clojure just kinda slow with large sequences? Reading some of the project Euler discussions people have calculated the first 10000 primes in other Lisps in under 100 miliseconds. I realize the JVM may slow things down and Clojure is relatively young, but I wouldn't expect a 100x difference.

Can someone enlighten me on a fast way to calculate prime numbers in Clojure?

16 Answers

Here is a Clojure solution. i is the current number being considered and p is a list of all prime numbers found so far. If division by some prime numbers has a remainder of zero, the number i is not a prime number and recursion occurs with the next number. Otherwise the prime number is added to p in the next recursion (as well as continuing with the next number).

(defn primes [i p]
  (if (some #(zero? (mod i %)) p)
    (recur (inc i) p)
    (cons i (lazy-seq (primes (inc i) (conj p i))))))
(time (do (doall (take 5001 (primes 2 []))) nil))
; Elapsed time: 2004.75587 msecs
(time (do (doall (take 10001 (primes 2 []))) nil))
; Elapsed time: 7700.675118 msecs

Update: Here is a much slicker solution based on this answer above. Basically the list of integers starting with two is filtered lazily. Filtering is performed by only accepting a number i if there is no prime number dividing the number with remainder of zero. All prime numbers are tried where the square of the prime number is less or equal to i. Note that primes is used recursively but Clojure manages to prevent endless recursion. Also note that the lazy sequence primes caches results (that's why the performance results are a bit counter intuitive at first sight).

(def primes
  (lazy-seq
    (filter (fn [i] (not-any? #(zero? (rem i %))
                              (take-while #(<= (* % %) i) primes)))
            (drop 2 (range)))))
(time (first (drop 10000 primes)))
; Elapsed time: 542.204211 msecs
(time (first (drop 20000 primes)))
; Elapsed time: 786.667644 msecs
(time (first (drop 40000 primes)))
; Elapsed time: 1780.15807 msecs
(time (first (drop 40000 primes)))
; Elapsed time: 8.415643 msecs

Idiomatic, and not too bad

(def primes
  (cons 1 (lazy-seq
            (filter (fn [i]
                      (not-any? (fn [p] (zero? (rem i p)))
                                (take-while #(<= % (Math/sqrt i))
                                            (rest primes))))
                    (drop 2 (range))))))
=> #'user/primes
(first (time (drop 10000 primes)))
"Elapsed time: 0.023135 msecs"
=> 104729

Plenty of answers already, but I have an alternative solution which generates an infinite sequence of primes. I was also interested on bechmarking a few solutions.

First some Java interop. for reference:

(defn prime-fn-1 [accuracy]
  (cons 2
    (for [i (range)
          :let [prime-candidate (-> i (* 2) (+ 3))]
          :when (.isProbablePrime (BigInteger/valueOf prime-candidate) accuracy)]
      prime-candidate)))

Benjamin @ https://stackoverflow.com/a/7625207/3731823 is primes-fn-2

nha @ https://stackoverflow.com/a/36432061/3731823 is primes-fn-3

My implementations is primes-fn-4:

(defn primes-fn-4 []
  (let [primes-with-duplicates
         (->> (for [i (range)] (-> i (* 2) (+ 5))) ; 5, 7, 9, 11, ...
              (reductions
                (fn [known-primes candidate]
                  (if (->> known-primes
                           (take-while #(<= (* % %) candidate))
                           (not-any?   #(-> candidate (mod %) zero?)))
                   (conj known-primes candidate)
                   known-primes))
                [3])     ; Our initial list of known odd primes
              (cons [2]) ; Put in the non-odd one
              (map (comp first rseq)))] ; O(1) lookup of the last element of the vec "known-primes"

    ; Ugh, ugly de-duplication :(
    (->> (map #(when (not= % %2) %) primes-with-duplicates (rest primes-with-duplicates))
         (remove nil?))))

Reported numbers (time in milliseconds to count first N primes) are the fastest from the run of 5, no JVM restarts between experiments so your mileage may vary:

                     1e6      3e6

(primes-fn-1  5)     808     2664
(primes-fn-1 10)     952     3198
(primes-fn-1 20)    1440     4742
(primes-fn-1 30)    1881     6030
(primes-fn-2)       1868     5922
(primes-fn-3)        489     1755  <-- WOW!
(primes-fn-4)       2024     8185 

If you don't need a lazy solution and you just want a sequence of primes below a certain limit, the straight forward implementation of the Sieve of Eratosthenes is pretty fast. Here's my version using transients:

(defn classic-sieve
  "Returns sequence of primes less than N"
  [n]
  (loop [nums (transient (vec (range n))) i 2]
    (cond
     (> (* i i) n) (remove nil? (nnext (persistent! nums)))
     (nums i) (recur (loop [nums nums j (* i i)]
                       (if (< j n)
                         (recur (assoc! nums j nil) (+ j i))
                         nums))
                     (inc i))
     :else (recur nums (inc i)))))
Related