エラトステネスの篩を使ってClojureで素数を求める

はじめに

Project Eulerの問題7, 素数10001番目を求める問題。
今回はエラトステネスの篩を使って実装してみた。

By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.
What is the 10 001st prime number?
Problem 7 - Project Euler

ソース

解説

関数sieveは篩という意味で、引数の数列の先頭の数値で、先頭以降の数列から先頭の数値の倍数になっているものをふるい落した数列を返却しています。この関数をprimesの中でiterateに渡して結果の先頭を取得すれば無限素数列となっているはずです。
ちなみに、primesが関数になっているのは、頭を抱え込まない*1ようにするためです。

試してみる

$ clj
Clojure 1.4.0
user=> (load-file "primes.clj")
#'user/primes
user=> (take 10 (primes))
(2 3 5 7 11 13 17 19 23 29)
user=> (take 100 (primes))
(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 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541)
user=> (nth (primes) 10000)
OutOfMemoryError Java heap space  clojure.core/filter (core.clj:2463)

どうもヒープを食い尽くすようです…アルゴリズムの検討が必要みたい。
また、改良して行ってみたいと思います。

ちなみに、clojure.contribにprimesがあるので、答えを調べてみる

$ clj
Clojure 1.4.0
user=> (use '[clojure.contrib.lazy-seqs :only (primes)])
nil
user=> (take 10 primes)
(2 3 5 7 11 13 17 19 23 29)
user=> (take 100 primes)
(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 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541)
user=> (nth primes 10000)
104743
user=>

ということで、104743らしいです。nthで10000を指定しているのは、インデックスが0番目から開始されるためです。

さらに別解

java.math.BigInteger に #nextProbablePrime() というメソッドがあって*2、これを使っても素数が求められそうです。

$ clj
Clojure 1.4.0
user=> (defn nextPrime [p] (.nextProbablePrime p))
#'user/nextPrime
user=> (def primes (iterate nextPrime (BigInteger. "2")))
#'user/primes
user=> (take 10 primes)
(2 3 5 7 11 13 17 19 23 29)
user=> (take 100 primes)
(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 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541)
user=> (nth primes 10000)
104743