在Clojure中快速生成质数

49
我一直在用Clojure解决Project Euler问题,以提高自己的能力,在这个过程中,我已经遇到了几次质数生成的问题。我的问题是这需要太长时间了。我希望有人能帮我找到一种高效的Clojure方法来解决这个问题。
起初我使用暴力算法,这很容易实现。但是,在Xeon 2.33GHz上计算10001个质数需要2分钟,这对规则来说太长了,而且总体上也太长了。以下是该算法:
(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))])))))

通过使用一个新的程序,它考虑了一些额外的规则(如6n +/- 1属性),取代next-prime-slow,我能够将速度提高到大约70秒。
接下来,我尝试使用纯Clojure制作Eratosthenes筛。我不认为我已经解决了所有的错误,但是我放弃了,因为它实在是太慢了(甚至比上面的还要糟糕)。
(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))))))))))

这很糟糕。如果数字150000更小,它也会导致堆栈溢出。尽管我使用了recur,但这可能是我的错。

接下来,我尝试了一个筛法,使用Java ArrayList上的Java方法。那花了相当长的时间和内存。

我最新的尝试是使用Clojure哈希映射的筛法,将所有数字插入筛子中,然后取消不是质数的数字。最后,它采用键列表,这些键是它找到的质数。查找10000个素数大约需要10-12秒。我不确定它是否完全调试好了。它也是递归的(使用recur和loop),因为我想成为Lispy。

所以对于这些问题,问题10(将2000000以下的所有质数相加)让我感到困扰。我的最快代码得出了正确的答案,但是它花费了105秒才能完成,并且需要相当多的内存(我给了它512 MB,只是为了不必担心它)。我的其他算法太慢了,我总是先停止它们。

我可以使用筛法在Java或C中快速计算出这么多质数,而且不需要使用太多内存。我知道我一定在我的Clojure / Lisp风格中漏掉了什么,导致出现问题。

我做错了什么吗?Clojure对于大序列是否有点慢?阅读一些项目Euler的讨论,人们在其他Lisps中计算出前10000个质数不到100毫秒。我知道JVM可能会减慢速度,而Clojure相对较年轻,但我不希望有100倍的差异。

有人能告诉我在Clojure中快速计算质数的方法吗?


你是想生成大量的质数,还是大质数?测试素性?目标是什么? - JP Alioto
我正在寻找一种通用算法。部分原因是为了提高我对语言的理解。一个问题要求找到第10001个质数,另一个问题要求计算小于2000000的所有数字之和。我预计还会有更多类似的问题。我上面的算法都是针对按顺序生成质数的。 - MBCook
不是答案,但我发现很有趣...http://bigdingus.com/2008/07/01/finding-primes-with-erlang-and-clojure/ - JP Alioto
我在Project Euler和Haskell中遇到了同样的问题,虽然不是同样的程度。我会在C和Haskell中实现相同的算法,而C程序只需要半秒钟,而Haskell则需要30秒。这主要是因为我不知道如何在Haskell中添加严格性,因为有些算法在两种语言中花费的时间大致相等。 - Dietrich Epp
这正是我遇到的问题。有些东西(涉及大量乘法、加法和其他操作)运行得非常快。但这个却非常慢。我正在努力学习是什么原因导致了这种情况,以便我可以修复它或避免它。 - MBCook
1
检查Alex Martelli的Python版本:https://dev59.com/snI95IYBdhLWcg3w-DH0 不同之处在于事先不知道会要求多少个数字。 - Hamish Grubijan
16个回答

29

这是另一种赞扬 Clojure的Java互操作性 的方法。在2.4 GHz Core 2 Duo上单线程运行,耗时374毫秒。我让Java中的BigInteger#isProbablePrime使用高效的Miller-Rabin实现来处理素性检查。

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

对于比这更大的数字,Miller-Rabin 的可靠性可能不是很高。该可靠性相当于 96.875% 的概率是质数(1 - .5^certainty)


2
哇,我甚至不知道BigInteger.isProbablePrime的存在。 - MBCook
4
1不是质数,BigInteger/isProbablePrime 在这一点上是正确的! - Omar Antolín-Camarena

22

我知道这是一个很旧的问题,但最近我也在寻找同样的东西,而这里的链接并不是我想要的(尽可能限制到函数类型,并且可以惰性生成我想要的每个质数)。

我偶然发现了一个很好的F# 实现,所以所有的功劳归他。我只是将它移植到 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)))

使用非常简单

(take 5 (gen-primes))    

这篇关于原始 F# 实现的链接文章非常棒!它教会了我很多有关 Clojure 中 maps 的知识。谢谢! - Kevin
3
Clojure不能像Scheme那样识别本地函数定义。使用defn定义的函数始终是全局的。对于像此示例中的相互递归本地函数,请使用letfn - Stuart Sierra
1
链接到 F# 文章已经不存在了。这是archive.org的链接 - A Sammich
2
我更喜欢这个作为 def。这是我的修改版本:(def primes "无限的质数懒序列。" ((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))关键点:
  • m 表示 map,n 表示 number
  • 内联了 reinsert
  • ((fn name [args] ...) args) 模式
  • 通过将 (step (inc n)) 可视化来说明递归情况的对称性
- yurrriq
@ASammich 我已经更新了答案中的链接。谢谢! - Matthew Murdoch
显示剩余4条评论

12

虽然我来得很晚,但是我会举一个使用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))

在我的2014款Macbook Pro(2.3GHz Core i7)上运行,我得到了以下结果:

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

1
不错,但是你如何从位集返回数字? - Phil Cooper
1
这是一个很好的答案,而且非常快。哦,要回到数字,你可以这样做(take-while #(not (= % -1)) (iterate #(.nextSetBit theBitSet (inc %)) 2)) - Matt Daley

11

看这里的最后一个例子: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)

4
虽然这种方法适用于小数字,但即使是这种方法也容易出现StackOverflowError错误 - 我在大约第2000个元素左右就遇到了此类错误。 - matt b
这实际上是特纳筛法(即试除法筛法),而不是埃拉托斯特尼筛法。它是二次的,而推迟的大约是~n^1.5(在n个产生的质数中)。 - Will Ness

4

(旁注:)链接博客中的Clojure代码相当于Haskell中的最优试除法primes = 2 : filter (\n -> all ((> 0).rem n) $ takeWhile ((<= n) . (^2)) primes ) [3..] - Will Ness

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

使用10M的堆栈大小,在一台2.1Gz的MacBook上,我可以在约33秒内得到第1001个质数。


我在寻找第100,000个质数时遇到了堆栈溢出。 - redcartel
每个质数都要被比它小的每个质数除,这比试除法要糟糕得多,试除法只检查到平方根。真正的埃拉托斯特尼筛法需要一个具有快速随机访问查找的数据结构。请参见http://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf。 - Phob

3

我刚开始学习Clojure,发现在欧拉计划中经常会遇到这样的问题。我写了一个很快的试除法素数算法,但是在每次运行除法时,随着数字变大速度变得非常慢。

于是我重新开始尝试,这次使用筛法:

(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

我对速度非常满意:它在一个运行于2009年MBP的REPL上运行。它主要很快是因为我完全避免了惯用的Clojure语法,而是像猴子一样循环。另外,我使用瞬态向量来处理筛子,所以速度也快了4倍,而没有完全采用不可变性。

编辑:在得到Will Ness的几个建议/错误修复后,现在运行速度快多了。


在你的sieve-step中,如果primes在索引i处已经为nil,则不应该为i调用clense-the-sieveclense应该调用自身,如(clense primes step (* step step))。对于i > 2sieve-step可以直接调用(clense primes (* 2 i) (* i i)),并将i增加2。 :) - Will Ness
哦,糟糕,我在代码中的第一个建议可能在某个时候被重写了太多次。 - SCdF
所以最后一个建议实际上使它运行得更慢了。我仍然不太熟悉Clojure优化,但我猜这是因为我很懒,每次都要进行(= 2 i)检查,而实际上我应该将整个过程拆分成唯一的代码来处理第一个调用。 - SCdF
可能的话,当你完全只处理奇数时,你可以至少获得2倍的加速:i=i+2; step=2*i。 - Will Ness
好的,现在我会放手不管了,它已经足够快了;-) - SCdF

2

2

这里提供了一种Clojure解决方案。 i 是当前正在考虑的数字,p 是迄今为止找到的所有质数列表。如果除以一些质数有余数为零,则数字i不是质数,并在下一个数字上进行递归。否则,在下一次递归中将质数添加到p(同时继续处理下一个数字)。

(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

更新: 这里有一个更加流畅的解决方案,基于上面的答案。 基本上,从2开始的整数列表被惰性地过滤。过滤是通过仅接受没有除以余数为零的质数来执行的数字i。所有质数都尝试,其中质数的平方小于或等于i。 请注意,primes被递归使用,但Clojure设法防止无限递归。还要注意,惰性序列primes缓存结果(这就是为什么性能结果乍一看有点反直觉的原因)。
(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

2
请在您的答案中添加一些解释,以便其他人可以从中学习。 - Nico Haase
2
经验性增长阶数,拜托了! - Will Ness
我的猜测是二次的,即5000个质数大约需要810毫秒,这个时间对吗? - Will Ness

1

根据威尔的评论,这是我对 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))))))))

比较的初始提交:

这是我将这个素数生成器从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

编辑:

性能比较,其中postponed-primes使用迄今为止看到的质数队列而不是对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)

(参考这里的答案:https://dev59.com/LnNA5IYBdhLWcg3wcNbF#7625207。你的代码与之相比如何?) - Will Ness
哦,我明白了。这是在递归调用生成器和队列之间做出的权衡。队列存储(最坏情况下)从 pp^2 之间的所有质数,但具有恒定时间的入队和出队操作。生成器的递归调用会生成从 2 到 c^.5 的质数,我想这很快就会减少。经重新考虑,我认为我更愿意采用递归调用的方法。 - Shaun
好的,经过进一步测试,递归调用和队列之间的性能差异似乎相当微小。对于较小的质数(可能是前500,000个质数),队列似乎稍微快一些,而随着质数数量的增加,递归调用可能会更占优势。或者也许我的机器在运行某些测试时同时在做其他事情。 - Shaun
我认为我已经将lazy-seq放在了正确的位置。例如,参见http://grimoire.arrdem.com/1.6.0/clojure.core/lazy_DASH_seq/。 - Shaun
同时,您关于将 q 作为 prime-help 的参数进行添加的提议是正确的。避免这些不必要的乘法运算可以使其运行速度快大约20%。 - Shaun
显示剩余9条评论

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接