#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 の次の素数 n を p の先頭に加えたもの」を自身に与えて
得られるリスト』の先頭に 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］
があるとすると, f に 2 と空のベクタ [] を与えると素数列が得られます.
ただし, 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 のソースコードからコピーしてきて 解答 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% ではないが, ある確率以上で素数であるか どうかを判定できるような判定法です.
primes-iv で頭に 2 を付け忘れてたので修正.