48

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?

Thumbnail
  • 12,695
  • 2
  • 24
  • 36
MBCook
  • 13,812
  • 7
  • 35
  • 41
  • Are you trying to generate lots of primes, large primes? Test primality? What's the goal? – JP Alioto Jun 07 '09 at 01:37
  • I was looking for a general algorithm. Partly this is just to improve my understanding of the language. One problem asked for the 10001st prime, one for the sum of all under 2000000. I expect there will be more. My algorithms above are all targeted at generating primes in order. – MBCook Jun 07 '09 at 01:43
  • Not an answer, but something I found interesting ... http://bigdingus.com/2008/07/01/finding-primes-with-erlang-and-clojure/ – JP Alioto Jun 07 '09 at 02:05
  • I had the same problem with Project Euler and Haskell, though not of the same magnitude. I'd implement the same algorithm in C and Haskell, and the C program would take a half second whereas Haskell took thirty. This is mostly due to the fact that I don't really know how to add strictness to Haskell, as some algorithms take about equal times in both languages. – Dietrich Epp Jun 07 '09 at 02:08
  • That's exactly what I'm running into. Some thing I've done (involving tons of multiplication and addition and other things) run blindingly fast. But this is horribly slow. I'm trying to learn what causes that so I can fix it or avoid it. – MBCook Jun 07 '09 at 02:14
  • 1
    Check Alex Martelli's Python version: http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python The difference is that one would not know how many numbers will be asked for in advance. – Hamish Grubijan Sep 20 '10 at 22:07

16 Answers16

29

Here's another approach that celebrates Clojure's Java interop. This takes 374ms on a 2.4 Ghz Core 2 Duo (running single-threaded). I let the efficient Miller-Rabin implementation in Java's BigInteger#isProbablePrime deal with the primality check.

(def certainty 5)

(defn prime? [n]
      (.isProbablePrime (BigInteger/valueOf n) certainty))

(concat [2] (take 10001 
   (filter prime? 
      (take-nth 2 
         (range 1 Integer/MAX_VALUE)))))

The Miller-Rabin certainty of 5 is probably not very good for numbers much larger than this. That certainty is equal to 96.875% certain it's prime (1 - .5^certainty)

Alexei - check Codidact
  • 17,850
  • 12
  • 118
  • 126
Ted Pennings
  • 1,001
  • 1
  • 14
  • 15
21

I realize this is a very old question, but I recently ended up looking for the same and the links here weren't what I'm looking for (restricted to functional types as much as possible, lazily generating ~every~ prime I want).

I stumbled upon a nice F# implementation, so all credits are his. I merely ported it to Clojure:

(defn gen-primes "Generates an infinite, lazy sequence of prime numbers"
  []
  (letfn [(reinsert [table x prime]
             (update-in table [(+ prime x)] conj prime))
          (primes-step [table d]
             (if-let [factors (get table d)]
               (recur (reduce #(reinsert %1 d %2) (dissoc table d) factors)
                      (inc d))
               (lazy-seq (cons d (primes-step (assoc table (* d d) (list d))
                                              (inc d))))))]
    (primes-step {} 2)))

Usage is simply

(take 5 (gen-primes))
Will Ness
  • 62,652
  • 8
  • 86
  • 167
Benjamin Podszun
  • 9,086
  • 3
  • 31
  • 41
  • The linked article about the original F# implementation is excellent! This taught me a lot about maps in Clojure. Thanks! – Kevin Mar 25 '15 at 19:51
  • 3
    Clojure does not recognize local function definitions like Scheme. `defn` is always global. For mutually-recursive local functions as in this example, use `letfn`. – Stuart Sierra May 07 '15 at 14:53
  • 1
    The linked to F# article is gone. Here's the [archive.org link](https://web.archive.org/web/20150710134640/http://diditwith.net/2009/01/20/YAPESProblemSevenPart2.aspx) – A Sammich Aug 21 '15 at 07:25
  • 2
    I prefer this as a `def`. Here's my tweaked version: `(def primes "Infinite, lazy sequence of prime numbers." ((fn step [m n] (or (some-> (get m n) (-> (->> (reduce #(update %1 (+ %2 n) conj %2) (dissoc m n))) (step (inc n)))) (-> (assoc m (* n n) (list n)) (step (inc n)) (->> (cons n) (lazy-seq))))) {} 2))` Key points: - `m` for map, `n` for number - Inlined `reinsert` - `((fn name [args] ...) args)` pattern - Illustrates the symmetry of the recursion cases by visually isolating `(step (inc n))` – yurrriq Aug 21 '15 at 09:47
  • @ASammich I've updated the link in the answer. Thanks! – Matthew Murdoch Aug 28 '15 at 06:53
  • Doesn't seem to terminate for me – RaymondMachira Jan 08 '16 at 16:12
  • @RaymondMachira It's a lazy, infinite sequence, so yeah, it won't terminate. You have to do something with the result, like `(take 5 l)` or `(nth l 50)` – Josh Oct 05 '16 at 18:43
  • This is VERY old by now, but still wanted to say thanks for the post as it helped me understand the algorithm better after having read the paper on it. I did my own version which is similar to yours and it helped my understand some things with lazy sequences that have proved invaluable. – Josh Oct 05 '16 at 18:55
  • `(update-in table [(+ prime x)] conj prime)` is just `(update table (+ prime x) conj prime)`. You could also use a [transient](https://clojure.org/reference/transients) map to speed it up a bit (untested). – Thumbnail Jul 27 '20 at 16:30
12

Very late to the party, but I'll throw in an example, using Java BitSets:

(defn sieve [n]
  "Returns a BitSet with bits set for each prime up to n"
  (let [bs (new java.util.BitSet n)]
    (.flip bs 2 n)
    (doseq [i (range 4 n 2)] (.clear bs i))
    (doseq [p (range 3 (Math/sqrt n))]
      (if (.get bs p)
        (doseq [q (range (* p p) n (* 2 p))] (.clear bs q))))
    bs))

Running this on a 2014 Macbook Pro (2.3GHz Core i7), I get:

user=> (time (do (sieve 1e6) nil))
"Elapsed time: 64.936 msecs"
Edward Dorrington
  • 121
  • 1
  • 2
11

See the last example here: http://clojuredocs.org/clojure_core/clojure.core/lazy-seq

;; An example combining lazy sequences with higher order functions
;; Generate prime numbers using Eratosthenes Sieve
;; See http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
;; Note that the starting set of sieved numbers should be
;; the set of integers starting with 2 i.e., (iterate inc 2) 
(defn sieve [s]
  (cons (first s)
        (lazy-seq (sieve (filter #(not= 0 (mod % (first s)))
                                 (rest s))))))

user=> (take 20 (sieve (iterate inc 2)))
(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71)
Will Ness
  • 62,652
  • 8
  • 86
  • 167
Wint
  • 2,128
  • 2
  • 20
  • 30
  • 4
    While this works for small numbers, even this method seems susceptible to StackOverflowError - I get them around the 2000th element or so – matt b Mar 10 '14 at 20:24
  • this is actually the *non-postponed sieve of Turner* (i.e. trial division sieve), not Eratosthenes. it's quadratic, while [the postponed](https://stackoverflow.com/a/8871918/849891) is about ~n^1.5 (in *n* primes produced) – Will Ness Feb 24 '20 at 15:22
4

Here's a nice and simple implementation:

http://clj-me.blogspot.com/2008/06/primes.html

... but it is written for some pre-1.0 version of Clojure. See lazy_seqs in Clojure Contrib for one that works with the current version of the language.

Will Ness
  • 62,652
  • 8
  • 86
  • 167
Jouni K. Seppänen
  • 36,462
  • 5
  • 69
  • 97
  • (side note:) the linked blog's clojure code amounts to `primes = 2 : filter (\n -> all ((> 0).rem n) $ takeWhile ((<= n) . (^2)) primes ) [3..]` in Haskell, the optimal trial division. – Will Ness Feb 27 '20 at 08:31
3

So I've just started with Clojure, and yeah, this comes up a lot on Project Euler doesn't it? I wrote a pretty fast trial division prime algorithm, but it doesn't really scale too far before each run of divisions becomes prohibitively slow.

So I started again, this time using the sieve method:

(defn clense
  "Walks through the sieve and nils out multiples of step"
  [primes step i]
  (if (<= i (count primes))
    (recur 
      (assoc! primes i nil)
      step
      (+ i step))
    primes))

(defn sieve-step
  "Only works if i is >= 3"
  [primes i]
  (if (< i (count primes))
    (recur
      (if (nil? (primes i)) primes (clense primes (* 2 i) (* i i)))
      (+ 2 i))
    primes))

(defn prime-sieve
  "Returns a lazy list of all primes smaller than x"
  [x]
  (drop 2 
    (filter (complement nil?)
    (persistent! (sieve-step 
      (clense (transient (vec (range x))) 2 4) 3)))))

Usage and speed:

user=> (time (do (prime-sieve 1E6) nil))
"Elapsed time: 930.881 msecs

I'm pretty happy with the speed: it's running out of a REPL running on a 2009 MBP. It's mostly fast because I completely eschew idiomatic Clojure and instead loop around like a monkey. It's also 4X faster because I'm using a transient vector to work on the sieve instead of staying completely immutable.

Edit: After a couple of suggestions / bug fixes from Will Ness it now runs a whole lot faster.

SCdF
  • 51,261
  • 23
  • 74
  • 108
  • in your `sieve-step`, you shouldn't call `clense-the-sieve` for `i` if `primes` at index `i` is already `nil`. `clense` should call itself as `(clense primes step (* step step))`. for `i > 2`, `sieve-step` can call directly `(clense primes (* 2 i) (* i i))`, and increment `i` by 2. :) – Will Ness Apr 29 '13 at 18:43
  • Oh whoops, I had that first suggestion in the code at some point, I guess too many rewrites.. – SCdF Apr 30 '13 at 01:22
  • So the last suggestion actually had it run slower. I'm still not really down with Clojure optimisation, but I'm guessing it's because I was lazy and had it do the (= 2 i) check each time, when really I should have split that whole thing off into unique code for that first call. – SCdF Apr 30 '13 at 01:57
  • probably. :) you can have at least 2x speedup when going fully to work on odds only: i=i+2; step=2*i. – Will Ness Apr 30 '13 at 01:59
  • OK, *now* I'm going to leave this alone, it's fast enough ;-) – SCdF Apr 30 '13 at 02:29
3
(defn sieve
  [[p & rst]]
  ;; make sure the stack size is sufficiently large!
  (lazy-seq (cons p (sieve (remove #(= 0 (mod % p)) rst)))))

(def primes (sieve (iterate inc 2)))

with a 10M stack size, I get the 1001th prime in ~ 33 seconds on a 2.1Gz macbook.

Joost
  • 31
  • 1
  • I got a stack overflow looking for the 100,000th prime – redcartel Jul 12 '12 at 20:50
  • Every prime is checked for division by every prime less than it. That's far worse than trial division, which only checks up to the square root. The true sieve of Eratosthenes requires a data structure with fast random access lookup. See http://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf – Phob Sep 22 '12 at 05:30
2

Here's a simple sieve in Scheme:

http://telegraphics.com.au/svn/puzzles/trunk/programming-in-scheme/primes-up-to.scm

Here's a run for primes up to 10,000:

#;1> (include "primes-up-to.scm")
; including primes-up-to.scm ...
#;2> ,t (primes-up-to 10000)
0.238s CPU time, 0.062s GC time (major), 180013 mutations, 130/4758 GCs (major/minor)
(2 3 5 7 11 13...
qu1j0t3
  • 650
  • 8
  • 12
2

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
wedesoft
  • 2,211
  • 23
  • 22
1

Based on Will's comment, here is my take on postponed-primes:

(defn postponed-primes-recursive
  ([]
     (concat (list 2 3 5 7)
             (lazy-seq (postponed-primes-recursive
                        {}
                        3
                        9
                        (rest (rest (postponed-primes-recursive)))
                        9))))
  ([D p q ps c]
     (letfn [(add-composites
               [D x s]
               (loop [a x]
                 (if (contains? D a)
                   (recur (+ a s))
                   (persistent! (assoc! (transient D) a s)))))]
       (loop [D D
              p p
              q q
              ps ps
              c c]
         (if (not (contains? D c))
           (if (< c q)
             (cons c (lazy-seq (postponed-primes-recursive D p q ps (+ 2 c))))
             (recur (add-composites D
                                    (+ c (* 2 p))
                                    (* 2 p))
                    (first ps)
                    (* (first ps) (first ps))
                    (rest ps)
                    (+ c 2)))
           (let [s (get D c)]
             (recur (add-composites
                     (persistent! (dissoc! (transient D) c))
                     (+ c s)
                     s)
                    p
                    q
                    ps
                    (+ c 2))))))))

Initial submission for comparison:

Here is my attempt to port this prime number generator from Python to Clojure. The below returns an infinite lazy sequence.

(defn primes
  []
  (letfn [(prime-help
            [foo bar]
            (loop [D foo
                   q bar]
              (if (nil? (get D q))
                (cons q (lazy-seq
                         (prime-help
                          (persistent! (assoc! (transient D) (* q q) (list q)))
                          (inc q))))
                (let [factors-of-q (get D q)
                      key-val (interleave
                               (map #(+ % q) factors-of-q)
                               (map #(cons % (get D (+ % q) (list)))
                                    factors-of-q))]
                  (recur (persistent!
                          (dissoc!
                           (apply assoc! (transient D) key-val)
                           q))
                         (inc q))))))]
    (prime-help {} 2)))

Usage:

user=> (first (primes))
2
user=> (second (primes))
3
user=> (nth (primes) 100)
547
user=> (take 5 (primes))
(2 3 5 7 11)
user=> (time (nth (primes) 10000))
"Elapsed time: 409.052221 msecs"
104743

edit:

Performance comparison, where postponed-primes uses a queue of primes seen so far rather than a recursive call to postponed-primes:

user=> (def counts (list 200000 400000 600000 800000))
#'user/counts
user=> (map #(time (nth (postponed-primes) %)) counts)
("Elapsed time: 1822.882 msecs"
 "Elapsed time: 3985.299 msecs"
 "Elapsed time: 6916.98 msecs"
 "Elapsed time: 8710.791 msecs"
2750161 5800139 8960467 12195263)
user=> (map #(time (nth (postponed-primes-recursive) %)) counts)
("Elapsed time: 1776.843 msecs"
 "Elapsed time: 3874.125 msecs"
 "Elapsed time: 6092.79 msecs"
 "Elapsed time: 8453.017 msecs"
2750161 5800139 8960467 12195263)
Community
  • 1
  • 1
Shaun
  • 81
  • 6
  • nice! I've edited to add the link. :) In Python the gain is much smaller; in your code the gain is very significant. Might be interesting to compare the [empirical orders of growth](http://en.wikipedia.org/wiki/Analysis_of_algorithms#Empirical_orders_of_growth) too, between the two versions, plain and postponed. – Will Ness Jul 13 '14 at 05:31
  • I'd pass `q` as one more argument to `prime-helper`, not to recalculate it all the time. Another thing to try is to allow multiplicity in the dictionary: { 45:[6,10] } kind of thing. Interesting, what might be the effect of that. The code will have to be more involved, but would that be worth it? – Will Ness Jul 13 '14 at 05:49
  • (cf. this answer here: http://stackoverflow.com/a/7625207/849891. how does your code fare in comparison?) – Will Ness Jul 13 '14 at 06:04
  • One difference between my postponed-primes and the one that you linked to is that, instead of using a "separate Primes Supply" `ps = (p for p in postponed_sieve())`, I'm initializing `ps` to a PersistentQueue with values 5 and 7, then pushing on primes as I see them with `(conj ps c)`. – Shaun Jul 13 '14 at 16:48
  • `user=> (time (nth (gen-primes) 10000)) "Elapsed time: 539.591 msecs" 104743` – Shaun Jul 14 '14 at 02:24
  • so your `postponed-primes` postpones adding the primes into dictionary (correctly) but it does needlessly add each prime to the tail of auxiliary list (queue). There should be a way to reuse a second invocation of the same function, `(postponed-primes)`, instead of using your PersistentQueue: `(prime-help {} 3 (postponed-primes) 9)`... (is this possible? advantageous? still interesting? :) :) ) – Will Ness Jul 14 '14 at 15:00
  • (also, if you could give times for the three functions also at 75,000, 50,000 and 25,000, it would enable us to compare them *algorithmically*) – Will Ness Jul 14 '14 at 15:03
  • I wrote it that way assuming that the second invocation of the `prime-help` function would needlessly force a recalculation of primes (and also assumed that the second call to `prime-help` would cause another call to `prime-help` and so on). Is there some sort of automatic memoization or other thing going on? – Shaun Jul 14 '14 at 15:43
  • I don't know Clojure to answer this. In Python, what works is a generator, maintaining its internal state, producing numbers one by one. First few primes (2,3,5,7) are specified inside. Each generator above, needs primes only up to the square root of its current production point, from the generator below. So the chain is finite and short enough. – Will Ness Jul 14 '14 at 17:37
  • Oh, I see. It's a tradeoff between the recursive call to the generator and the queue. The queue stores (worst case) all the primes between `p` and `p^2` but has constant time enqueuing and dequeuing operations. The recursive calls to the generator generate primes from 2 up to `c^.5`, which I guess tapers off quite quickly. After reconsideration, I think I'd rather take the recursive call approach. – Shaun Jul 15 '14 at 01:06
  • Ok, so from further testing, the performance difference between the recursive call and the queue seems to be pretty marginal. Queue seems to be a little faster for smaller numbers of primes (maybe the first 500,000 primes) while the recursive call probably gets the edge as you go past that. Or maybe my machine happened to be doing something else while running some of those tests. – Shaun Jul 16 '14 at 06:02
  • are you sure it should be `(if (< c q) (cons c (lazy-seq (prime-help D p ...` and not `(if (< c q) (lazy-seq (cons c (prime-help D p ...`? http://stackoverflow.com/a/7625207/849891 seems to be doing the latter. Maybe that's why its `gen-primes` is a little bit faster than your `primes`. -- With the new data we see that the empirical orders of growth for both primes and gen-primes are around n^1.25 (at the tested range), for postponed with queue - n^1.11..1.16 and with recursive call - n^1.05..1.15. Of course the numbers are approximate, but the general picture can be seen there. – Will Ness Jul 16 '14 at 08:49
  • I think I have the `lazy-seq` in the right place. See, e.g., http://grimoire.arrdem.com/1.6.0/clojure.core/lazy_DASH_seq/. – Shaun Jul 17 '14 at 05:32
  • Also, you were right about adding `q` as an argument to prime-help. Avoiding those unnecessary multiplications makes it about 20% faster. – Shaun Jul 17 '14 at 05:33
1

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
Thumbnail
  • 12,695
  • 2
  • 24
  • 36
Tom
  • 21
  • 1
  • 1
0

From: http://steloflute.tistory.com/entry/Clojure-%ED%94%84%EB%A1%9C%EA%B7%B8%EB%9E%A8-%EC%B5%9C%EC%A0%81%ED%99%94

Using Java array

(defmacro loopwhile [init-symbol init whilep step & body]
  `(loop [~init-symbol ~init]
     (when ~whilep ~@body (recur (+ ~init-symbol ~step)))))

(defn primesUnderb [limit]
  (let [p (boolean-array limit true)]
    (loopwhile i 2 (< i (Math/sqrt limit)) 1
               (when (aget p i)
                 (loopwhile j (* i 2) (< j limit) i (aset p j false))))
    (filter #(aget p %) (range 2 limit))))

Usage and speed:

user=> (time (def p (primesUnderb 1e6)))
"Elapsed time: 104.065891 msecs"
KIM Taegyoon
  • 1,638
  • 20
  • 14
0

I just started using Clojure so I don't know if it's good but here is my solution:

(defn divides? [x i]
  (zero? (mod x i)))

(defn factors [x]
    (flatten (map #(list % (/ x %)) 
                 (filter #(divides? x %) 
                        (range 1 (inc (Math/floor (Math/sqrt x))))))))

(defn prime? [x]
  (empty? (filter #(and divides? (not= x %) (not= 1 %)) 
                 (factors x))))

(def primes 
  (filter prime? (range 2 java.lang.Integer/MAX_VALUE)))

(defn sum-of-primes-below [n]
  (reduce + (take-while #(< % n) primes)))
Will Ness
  • 62,652
  • 8
  • 86
  • 167
user2428814
  • 1
0

After coming to this thread and searching for a faster alternative to those already here, I am surprised nobody linked to the following article by Christophe Grand :

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

As well as a lazy version :

(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)))))
nha
  • 16,039
  • 11
  • 76
  • 108
  • 1
    `(enqueue sieve candidate candidate)` can be safely replaced with `(enqueue sieve (- (* candidate candidate) candidate candidate) candidate)`. this will entail [change](https://stackoverflow.com/a/10733621/849891) that will bring memory requirements down from O(n) to O(sqrt(n)), and depending on how good the dictionary implementation in Clojure is, could even give a sizeable constant factor (or even slight algorithmic (i.e. empirical orders of growth) improvement). coincidentally, the linked answer links to 2002(?) Python recipe which is an exact equivalent of this one in your answer. :) – Will Ness Nov 04 '20 at 15:35
  • 1
    (see also: https://stackoverflow.com/search?q=user%3A849891+postponed+sieve) – Will Ness Nov 04 '20 at 15:36
0

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
NikoNyrh
  • 3,122
  • 1
  • 13
  • 30
0

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)))))
miner49r
  • 1,047
  • 1
  • 12
  • 10