Skip to content

Instantly share code, notes, and snippets.

@kohyama kohyama/README.md
Last active Oct 13, 2015

Embed
What would you like to do?
Project Euler Problem 10

Project Euler Problem 10 解答とコメント

#clojure 入門者向け勉強会 #mitori_clj 第二週担当分

Project Euler Problem 10
2,000,000 未満の素数の総和を求めよ. とのことです.

結果を出すために手に入るリソースを有効活用するという考えのもとでは --この問題自体の解答を述べている 他の方のコードを利用させていただくことを除けば--, clojure 1.2 時代に clojure.contrib に存在していた素数列 primes の定義を利用するのが正解と思います.

primes と, その定義に必要な defvar の定義をそのままコピーしてきて, (apply + (take-while #(< % 2000000) primes)) で 範囲内の和を求めるコードを解答として pep010.clj に記述しておきます.

中身のコードについては, 完全には理解していませんが, 後で分かる範囲で解説します.

ところでこの primes を含む lazy-seqs.clj ですが clojure 1.3 以後, どこに行ったのか分かりません. 誰か知っていたら教えてください.

さて, 今回の目的は答えを出す事ではなく, この問題を題材に勉強することですので, 自分で少しやってみます.

素数列が得られれば, 和を求める部分は (apply + (take-while #(< % 2000000) 素数列)) で問題ないと思いますので, 素数列を求めるところにフォーカスを当てます.

記述量が少ないところから始めたいと思います.

無限数列に対して篩をかける

『先頭の値が h, 以後が r であるような数列を与えられると, 「数列 r から, 値 h で割切れる数を除いた数列を自身に与えて得られる数列」 の先頭に値 h を付加した数列を返す手続き f』があるとすると, f に, 2 から始まり 1ずつ増加する自然数列を与えて得られるものが素数列です.

  • 必要になるまで無限数列の残りを計算しないように cons の第二引数を遅延シーケンスにすること
  • シーケンス全体への参照を保持しないよう関数として記述すること,

に注意して再帰で記述します.

(defn gen-primes-i [] 
  ((fn f [[h & r]]
     (cons h (lazy-seq (f (remove #(zero? (rem % h)) r)))))
   (iterate inc 2)))

100 未満の素数を求めてみます.

user=> (take-while #(< % 100) (gen-primes-i))
(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97)

問題文中の「10未満の素数の和」で動作確認します.

user=> (apply + (take-while #(< % 10) (gen-primes-i)))
17

いいですね. では, 2,000,000 未満の素数の和を求めてみます.

user=> (apply + (take-while #(< % 2000000) (gen-primes-i)))
StackOverflowError   clojure.lang.LazySeq.sval (LazySeq.java:42)

残念. スタックが溢れてしまいます.

環境にもよると思いますが, 私の手元では 10,000 未満の素数の総和は求まりますが, 100,000 未満の素数の総和は求まりません.

user=> (apply + (take-while #(< % 10000) (gen-primes-i)))
5736396
user=> (apply + (take-while #(< % 100000) (gen-primes-i)))
StackOverflowError   clojure.core/seq (core.clj:133)

もう少しがんばってみましょう.

既知の素数のリストから次の素数を順に求める

2 から h までの既知の素数が大きい順に入ったリスト p を受け取り, 『「h の次の素数 np の先頭に加えたもの」を自身に与えて 得られるリスト』の先頭に h を追加する手続き f]があるとすると, f に要素 2 のみからなるリストを与えると素数列が得られます.

ただし, h の次の素数 n は, h + 1 から始まり, 1ずつ増加する自然数列の要素 c を候補として順に見て, p の中のいずれかの要素 % で, c が割り切れる間は, 素数では無いとして捨てていって, 最初に現れた数です.

(require '[clojure.test :refer :all])

(defn gen-primes-ii []
  ((fn f [[h & _ :as p]]
     (let [n (first (drop-while (fn [c] (some #(zero? (rem c %)) p))
                                (iterate inc (inc h))))]
       (cons h (lazy-seq (f (cons n p))))))
   '(2)))

(is (= (take-while #(< % 100) (gen-primes-ii))
      '(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73
        79 83 89 97)))

(is (= (apply + (take-while #(< % 10) (gen-primes-ii))) 17))

2,000,000未満の素数の総和を求めてみます.

user=> (apply + (take-while #(< % 2000000) (gen-primes-ii)))

手元の環境では 5分以上放置して答えが出ないようなので中断. もう少し効率化してみます.

平方根以下の素数で試し割り

ある自然数 c が素数であるかどうかは, c 未満の素数全てで割り切れないことを確認する必要はなく, c の正の平方根以下の素数 全てで割り切れないことを確認すれば十分です.

試し割りする素数をある数以下に限定して効率を上げたいので, それまでに求まった素数を小さい順で保持するようにします.

[素数 m と, m 未満の既知の素数が小さい順に入ったベクタ p を受け取り, 『「m の次の素数 n」と「p の末尾に m を加えたベクタ q」 を自身に与えて得られるリスト』の先頭に m を追加する手続き f] があるとすると, f2 と空のベクタ [] を与えると素数列が得られます.

ただし, m の次の素数 n は, m + 1 から始まり, 1ずつ増加する自然数列の要素 c を候補として順に見て, q の要素 % が, c の正の平方根以下である数列 の中のいずれかの要素 % (先の % とは違うので注意) で, c が割り切れる間は, 素数では無いとして捨てていって, 最初に現れた数です.

(require '[clojure.test :refer :all])

(defn gen-primes-iii []
  ((fn f [m p]
     (let [q (conj p m)
           n (first
               (drop-while
                 (fn [c] (some #(zero? (rem c %)) 
                               (take-while #(<= (* % %) c) q)))
                 (iterate inc (inc m))))]
       (cons m (lazy-seq (f n q)))))
   2 []))

(is (= (take-while #(< % 100) (gen-primes-iii))
      '(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73
        79 83 89 97)))

(is (= (apply + (take-while #(< % 10) (gen-primes-iii))) 17))

(apply + (take-while #(< % 2000000) (gen-primes-iii)))

手元の環境では計算時間は約 10 秒でした.

user=> (time (apply + (take-while #(< % 2000000) (gen-primes-iii))))
"Elapsed time: 10178.917 msecs"

3 以上の素数は全て奇数であることを利用し,

(defn gen-primes-iv []
  (cons 2
    ((fn f [m p]
       (let [q (conj p m)
             n (first
                 (drop-while
                   (fn [c] (some #(zero? (rem c %))
                                 (take-while #(<= (* % %) c) q)))
                   (drop 1 (iterate #(+ % 2) m))))]
         (cons m (lazy-seq (f n q)))))
     3 [2])))

とすると, ちょっとだけ速くなります.

user=> (time (apply + (take-while #(< % 2000000) (gen-primes-iv))))
"Elapsed time: 8861.406 msecs"

clojure.contrib の primes の解説

clojure.contrib のソースコードからコピーしてきて 解答 pep010.clj に掲載した primes について分かる範囲で説明します.

Recursion on the cycle of gaps の理論を利用しているようです.

[2 3 5 7] までは既知として G(7#) を利用しています.

wheel (cycle [2 4 2 4 6 2 6 4 2 4 6 6 2 6  4  2
              6 4 6 8 4 2 4 2 4 8 6 4 6 2  4  6
              2 6 6 4 2 4 6 2 6 4 2 4 2 10 2 10])

の部分が G(7#) の第二要素から始めて循環させたものです.

上記ページには, G(pk#) を順に更新すれば, G(pk-1#) の最初の要素に 1 を足したものが次の素数 pk になると書いてありますが, primes のコードでは, G(pk#) を更新することはせず, G(7#) を使って, 次の素数候補を絞り,

(some #(zero? (rem n %))
      (take-while #(<= (* % %) n) primes))

と試し割りをしているようです.

ただ, G(pk) が以後の素数の候補として どう使えるのかについては, 上記ページからは読み取れず, 理論との対応がよく分かりませんでした.

どなたかお分かりになったら教えてください.

(some #(zero? (rem n %))
      (take-while #(<= (* % %) n) primes))

の箇所で, もうひとつ面白いのは, 試し割りするための既知の素数列を覚えおく機構を用意せず, 自分自身を呼び出しているところです.

かっこいいですが, 下手にまねをすると, 無限ループに陥ります.

というか陥りました.

うまくまねできなかったので, gen-primes-iii では, 求めた素数を追加しながら既知の素数を順に引数で渡しています.

あと, defvar は doc-string やメタデータの追加を書きやすくするマクロのようです.

ちなみに計算時間は, 手元の環境で

user=> (time (apply + (take-while #(< % 2000000) primes)))
"Elapsed time: 11260.111 msecs"

となり, gen-primes-iii (ともちろん gen-primes-iv) の方が少しだけ速いです.

その他の方法

試し割りはしない

無限数列に篩をかける方法で, 試しに割って 0 になる数を除くのではなく, 倍数列を除くのが本来の「エラストテネスの篩」であるとし, 無限数列から無限数列を除く方法でこれを実現した 「エラストテネスの無限の篩」は, clojure らしくてとても素敵です.

素数判定ありき

素数であるかどうかという判定は, 公開鍵暗号システムなどでとても実用的に重要な役割を果たします. Wikipedia/素数判定

素数であるか判定するアルゴリズムが先にある場合, 自然数列から判定に合格する数のみ抜き出すことによっても素数列は求まります.

ただし, 判定アルゴリズムが, その数の平方根以下の素数を利用する --または, 素数であるかに関わらず, 2以上平方根以下の全ての整数で割り切れないことを確認する-- ならば, 本質的には上記 gen-primes-iii と同じです.

RSA暗号で重要な役割を果たす非常に大きな素数の生成には, 確率的素数判定法である 「ミラー・ラビン法 (Wikipedia/ミラー-ラビン素数判定法)」 などが用いられるようです.

確率的素数判定法とは, 100% ではないが, ある確率以上で素数であるか どうかを判定できるような判定法です.

(require '[clojure.test :refer :all])
;; just copied from https://github.com/richhickey/clojure-contrib/blob/master/src/main/clojure/clojure/contrib/def.clj#L23
(defmacro defvar
"Defines a var with an optional intializer and doc string"
([name]
(list `def name))
([name init]
(list `def name init))
([name init doc]
(list `def (with-meta name (assoc (meta name) :doc doc)) init)))
;; just copied from https://github.com/richhickey/clojure-contrib/blob/master/src/main/clojure/clojure/contrib/lazy_seqs.clj#L66
(defvar primes
(concat
[2 3 5 7]
(lazy-seq
(let [primes-from
(fn primes-from [n [f & r]]
(if (some #(zero? (rem n %))
(take-while #(<= (* % %) n) primes))
(recur (+ n f) r)
(lazy-seq (cons n (primes-from (+ n f) r)))))
wheel (cycle [2 4 2 4 6 2 6 4 2 4 6 6 2 6 4 2
6 4 6 8 4 2 4 2 4 8 6 4 6 2 4 6
2 6 6 4 2 4 6 2 6 4 2 4 2 10 2 10])]
(primes-from 11 wheel))))
"Lazy sequence of all the prime numbers.")
(is (= (apply + (take-while #(< % 10) primes)) 17))
(apply + (take-while #(< % 2000000) primes))
@kohyama

This comment has been minimized.

Copy link
Owner Author

commented Dec 10, 2012

primes-iv で頭に 2 を付け忘れてたので修正.

@kohyama

This comment has been minimized.

Copy link
Owner Author

commented Dec 11, 2012

defvar について誤解していたようなので修正.

@kohyama

This comment has been minimized.

Copy link
Owner Author

commented Dec 11, 2012

関数版は gen- をプレフィクスするべきと思ったので修正.

@tnoda

This comment has been minimized.

Copy link

commented Dec 13, 2012

下から読んでます.

RSA暗号で重要な役割を果たす非常に大きな素数の生成には, 確率的素数判定法である 「ミラー・ラビン法(Wikipedia/ミラー-ラビン素数判定法)」 などが用いられるようです.

素の Clojure でもミラーラビン使えるんですよね.

(def prime? #(-> % str BigInteger. (.isProbablePrime 40)))

手元の環境で試したところ,この prime? 関数は 100 万未満のすべての素数を正しく判定できました.


以下余談.

OpenJDK7 の .isProbablePrime ではミラーラビンとルーカスレーマーとを併用していました.

http://www.docjar.com/html/api/java/math/BigInteger.java.html

Clojure を使うと「今ならミラーラビンだけでなくルーカスレーマーもついてきます」と聞くと何だかおトクな気分になれますね.

@ypsilon-takai

This comment has been minimized.

Copy link

commented Dec 14, 2012

.isProbablePrime これすごい。

@emanon001

This comment has been minimized.

Copy link

commented Dec 18, 2012

数学の苦手な僕でも分かりやすい解説で、大変勉強になりました。
他の問題でも素数列を使用するものが多いので、今の段階で知識を整理できて良かったです。
そして、エラストテネスの無限の篩、とてもよい。


.isProbablePrime
なんと、Java の資産恐るべし。

@kohyama

This comment has been minimized.

Copy link
Owner Author

commented Dec 19, 2012

.isProbablePrime 👍

@plaster

This comment has been minimized.

Copy link

commented Dec 19, 2012

私も数学苦手ですが、ギャップで次の候補を求めるの興味深いです。
リンク先を見ると G(7#) は 10 から始まっていますが、pep010.clj の wheel は 2 から始まっているのもよくわかっていません。もちろんこれでいきなり10だと、次が21なので明らかに変なのですが……先頭の10は、一番後ろに回っている?

「次の素数」(おそらく確定?)との関係も全然わかっていないのですが、素数が無限にあることを背理法で証明するときに「既知の素数すべての積+1」を持ち出したのを思い出しました。的外れだったらすみません。

そしてJavaすごいですね……isProbablePrime……

@kohyama

This comment has been minimized.

Copy link
Owner Author

commented Dec 19, 2012

ギャップ理論と関係あるのかないのか「3 より大きい全ての素数は必ず 6の倍数 - 1 か 6の倍数 + 1 であるので 2, 3, 6 - 1, 6 + 1, 2_6 - 1, 2_6 + 1, 3_6 - 1, 3_6 + 1,... を素数の候補として使う」というのも, どこかで見た事があります.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.