代码之家 › 专栏 › 技术社区 › MBCook

Clojure中的快速素数生成

primes lisp clojure

MBCook · 技术社区 · 15 年前

我一直在努力解决 Project Euler Clojure中的问题会变得更好,我已经遇到过几次质数生成。我的问题是时间太长了。我希望有人能帮我找到一个有效的方法来做这件事。

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

通过将next prime slow替换为考虑了一些额外规则(比如6n++1属性)的新例程,我可以将速度提高到大约70秒。

下一步我试着用纯正的Clojure做一个筛子。我不认为我把所有的虫子都弄出来了,但我放弃了,因为它太慢了(我想比上面更糟)。

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

接下来,我尝试了一个筛子,在Java ArrayList上使用Java方法。那花了不少时间和记忆。

我最新的尝试是使用Clojure散列映射筛选,插入筛选中的所有数字,然后分解不是素数的数字。最后,它得到了密钥列表,这是它找到的素数。找到10000个质数大约需要10-12秒。我不确定它是否已经完全调试好了。它也是递归的(使用recur和loop),因为我想保持口齿不清。

所以,有了这些问题,问题10(总结所有200万以下的素数)就要把我累死了。我最快的代码给出了正确的答案,但它花了105秒,并且需要相当多的内存(我给它512MB,只是为了不必费心)。我的其他算法花了很长时间,我总是先停止它们。

我可以用筛子快速计算Java或C语言中的许多素数,而且不需要占用太多内存。我知道我的Clojure/Lisp风格中一定缺少导致问题的东西。

我真的做错什么了吗?Clojure是不是太慢了?阅读一些关于欧拉计划的讨论,人们已经在100毫秒内计算出了其他Lisps中的前10000个素数。我意识到JVM可能会减慢速度,Clojure还比较年轻,但我预计不会有100倍的差异。

有人能教我一种快速计算Clojure中素数的方法吗?

15 回复 | 直到 5 年前

Alexei - check Codidact 8 年前

下面是另一种庆祝 Clojure's Java interop . 这需要374ms的2.4ghz核心2双工(运行单线程)。我让有效率的 Miller-Rabin Java的实现 BigInteger#isProbablePrime 处理素性检查。

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

这个对于比这个大得多的数字来说,5的确定性可能不是很好。这种确定性等于 96.875% 1 - .5^certainty )

Will Ness Derri Leahy 4 年前

我偶然发现一个 F# implementation ,所以所有的学分都是他的。我只是把它转给了克洛朱尔:

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

(take 5 (gen-primes))

Edward Dorrington 10 年前

派对很晚了,但我将提供一个使用Java位集的示例:

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

在2014 Macbook Pro(2.3GHz Core i7)上运行此程序,我得到:

user=> (time (do (sieve 1e6) nil))
"Elapsed time: 64.936 msecs"

Will Ness Derri Leahy 4 年前

请参见下面的最后一个示例: 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 Derri Leahy 8 年前

下面是一个简单的实现:

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

... 但它是为Clojure的1.0以前的版本编写的。见 lazy_seqs

Joost 13 年前

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

在10米的堆栈大小下,我在2.1Gz的macbook上用大约33秒的时间就能得到1001th prime。

SCdF 11 年前

所以我刚从Clojure开始,是的,这在Euler项目上出现了很多,不是吗?我写了一个相当快速的试算除法素数算法,但是在每次除法都变得非常缓慢之前,它并没有太大的扩展性。

所以我又开始了,这次用筛分法:

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

使用和速度:

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

我对这个速度很满意:它的REPL快用完了,运行速度是2009年的MBP。最快的是,因为我完全避开了惯用的Clojure,而是像猴子一样绕来绕去。它的速度也快了4倍,因为我使用了一个瞬态向量来处理筛子,而不是保持完全不变。

编辑:

qu1j0t3 14 年前

下面是一个简单的筛选方案:

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

这里有一个10000以下的质数:

#;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...

Community vonPryz 7 年前

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

首次提交供比较:

this prime number generator 从Python到Clojure。下面返回一个无限延迟序列。

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

用法:

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

编辑:

性能比较,其中使用到目前为止看到的素数队列,而不是递归调用 :

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)

KIM Taegyoon 11 年前

发件人: 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

使用Java数组

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

使用和速度:

user=> (time (def p (primesUnderb 1e6)))
"Elapsed time: 104.065891 msecs"

user2428814 9 年前

我刚开始使用Clojure,所以我不知道它是否好,但这里是我的解决方案:

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

nha 8 年前

在来到这个主题并寻找一个比这里已有的更快的替代方案之后,我很惊讶没有人与下面的链接 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))))))

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

NikoNyrh 6 年前

已经有很多答案了,但是我有一个替代的解决方案,它可以生成一个无穷的素数序列。我也有兴趣提出一些解决方案。

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

本杰明@ https://stackoverflow.com/a/7625207/3731823 是 primes-fn-2

国家卫生局@ https://stackoverflow.com/a/36432061/3731823 是 primes-fn-3

我的实现是 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?))))

报告的数字(计算前N个素数的时间,以毫秒为单位)是运行5次之后最快的,两次实验之间没有JVM重新启动,因此您的里程可能会有所不同:

                     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

miner49r 5 年前

如果你不需要一个懒惰的解决方案,你只需要一个素数序列低于某个限制,那么直接实现埃拉托舍内斯筛是相当快的。这是我用瞬变的版本:

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

Thumbnail 5 年前

很地道,也不错

(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