すごい暇だったのでエラトステネスの篩(Sieve of Eratosthenes)分かる範囲で調べてみようと思った。
数学において、エラトステネスの篩（エラトステネスのふるい）は、指定された整数以下の全ての素数を発見するための単純なアルゴリズムである。古代ギリシアの科学者、エラトステネスが考案したとされるため、この名がある。(エラトステネスの篩 - Wikipedia)
アルゴリズムの詳細についてはWikipediaを見てもらうか他の人の説明を検索してもらった方が良いと思うけど一応簡単におさらい。
- E1(リストの用意) 2以上N以下(Nは目的とする範囲の上限)の整数をリストに追加。
- E2(フラグの用意) リスト中全ての整数について「素数であるフラグ」を立てる。
- E3(リストの走査) リストの値が小さいものから順にフラグを調べる
素直に実装する
上述のアルゴリズムをそのままRで実装するとこんな感じになる。
### version 1 ## 普通に定義通りに sieve1 <- function(limit = 1e6, return = FALSE){ sieve <- logical(limit) sieve[1] <- TRUE # 非素数はTRUEで表す for(i in 2:sqrt(limit)){ # 上限の平方根まで探す if(!sieve[i]){ # 素数なら sieve[seq(i*2, limit, by=i)] <- TRUE # その倍数を取り除く } } prime <- (1:limit)[!sieve] if(return) prime }
logical()は初期値がFALSEなのでFALSEをフラグONとみなしている以外はほとんどそのまんま書いてあるようなもの。sieve1(100, TRUE)とすれば100以下の素数がベクトルで返ってくる。
1000万以下を探索範囲とすると大体3秒くらいで返ってくる(OS:Mac OSX, プロセッサ: 2.26 GHz Intel Core 2 Duo, メモリ: 4 GB 1067 MHz DDR3)。
そもそもエラトステネスの篩は早いので、普通に使う程度の素数(?)なら今みたいに適当に書いたって十分使える。
ちょっと工夫する
そのまんま実装するとちょっと無駄がある。例えば2を除けば素数は奇数なのだから、初めから奇数のみ調べるようにしたらいい。また、小さい方の値から調べていっているのだから、素数pの整数倍のうち、p*pより小さなものは既にフラグがオフになっている。フラグをオフにする際のループはp*pよりスタートすればよい。
というわけでちょっと無駄を省いて実装すると、こんな感じになる。
### version 2 ## やや改良 sieve2 <- function(limit = 1e6, return = FALSE){ sieve <- logical(limit) sieve[1] <- TRUE sieve[seq(4, limit, by=2)] <- TRUE # 4以上の偶数を取り除く for(i in seq(3, sqrt(limit), by=2 )){ # 奇数だけ探す if(!sieve[i]){ ## i^2より小さな非素数はiより小さい約数を持つので ## すでにリストから外れている => 探索はi^2から ## i^2は奇数、i^2+iは偶数、i^2+2iは奇数... ## => 偶数は既に外れているのでステップは2iずつで良い sieve[seq(i^2, limit, by=i*2)] <- TRUE } } prime <- (1:limit)[!sieve] if(return) prime }
1000万までの素数をおおよそ2秒で求められる。
もうちょっと工夫する
さっきは2の倍数を調べないという方策をとったが、それならいっそ2の倍数のフラグを用意しなければいい。そうすればメモリも有効活用できる。添字と値の関係に注意して実装するとこんな感じ。
### version 3 ## 2の倍数を最初から除く sieve3 <- function(limit = 1e6, return = FALSE){ ## 2i + 1 = nとし、3から開始する sieve <- logical((limit-1)/2) # 1を除く数の半分が上限 for(i in 1:(floor(sqrt(limit)-1)/2)){ if(!sieve[i]){ ## 数nのindexは(n-1)/2で求められる ## n = 2i + 1 なので n^2 = 4i^2 + 4i + 1, そのindexは ## ((4i^2 + 4i + 1) - 1)/2 = 2i(i + 1) ## 数xのindexは(x-1)/2 ## 数xに2nを足すと、x + 2n = x + 4i + 2、そのindexは ## ((x + 4i + 2) - 1)/2 = (x -1 + 4i + 2)/2 で2i+1増えてる ## よって2nを足すごとにindexは2i + 1ずつ増える sieve[seq(2*i*(i+1), length(sieve), by = 2*i+1)] <- TRUE } } prime <- c(2, (1:length(sieve))[!sieve] * 2 + 1) if(return) prime }
1000万までの素数は1.3秒くらいで求められる。
もっともっと工夫したい
3も素数だし、3の倍数も初めから除いてしまえばいい。そうすると調べるべきはまたは
の数のみになる。だいぶメモリが節約できる。
が、Rだと1オリジンだからそのへんに落ちてるソースそのまま使えないし、分岐や条件式が増えるとすぐ重くなるしでどうも今ひとつ。あとseq()がエラー吐くの防ぐためにちょっとその場しのぎでアレなやり方してる。
まあとにかくこんな感じになる。
### version 4 ## 2の倍数と3の倍数を最初から除く sieve4 <- function(limit = 1e6, return =FALSE){ sieve1 <- logical(limit/5) sieve1[1] <- TRUE sieve5 <- logical(limit/5) slen <- length(sieve1) slim <- ceiling(sqrt(limit)/6) for(i in 1:slim){ if(!sieve1[i]){ p <- 6 * i - 5 sieve1[seq(6*i^2 - 10*i + 5, slen, by = p)] <- TRUE sieve5[seq(6*i^2 - 6*i + 1, slen, by = p)] <- TRUE } if(!sieve5[i]){ p <- 6 * i - 1 sieve1[seq(6*i^2 - 2*i + 1, slen, by = p)] <- TRUE sieve5[seq(6 * i^2, slen, by = p)] <- TRUE } } prime <- sort(c(2, 3, seq(1, limit, by = 6)[!sieve1], seq(5, limit, by = 6)[!sieve5])) if(return) prime }
計算速度としてはversion3とほとんど変わらない。範囲が狭いとむしろ遅いくらい。
もっとなんか別の使おう
アトキンの篩(Sieve of Atkin)というのがあって、エラトステネスの篩より高速に動作するらしい。というわけで説明(Spaghetti Source - アトキンのふるい)を読みながら丸写ししてみたけどさっぱり速度がでない。多分条件判定減らすとかしてRに合った実装しないとダメなんだと思う。といって元論文読む能力がないのでmoudameda。
とりあえずこんな感じ。まあそりゃ遅いよなって感じがひしひしと。
### アトキンの篩 atkin <- function(limit = 1e6, return = FALSE){ isprime <- logical(limit) z <- 1 while(z <= 5){ y <- z while(y <= sqrt(limit)){ x <- 1 while(x <= sqrt(limit) && (n <- 4*x*x+y*y) <= limit){ isprime[n] <- !isprime[n] x <- x + 1 } x <- y + 1 while(x <= sqrt(limit) && (n <- 3*x*x-y*y) <= limit){ isprime[n] <- !isprime[n] x <- x + 2 } y <- y + 6 } z <- z + 4 } z <- 2 while(z <= 4){ y <- z while(y <= sqrt(limit)){ x <- 1 while(x <= sqrt(limit) && (n <- 3*x*x+y*y) <= limit){ isprime[n] <- !isprime[n] x <- x + 2 } x <- y + 1 while(x <= sqrt(limit) && (n <- 3*x*x-y*y) <= limit){ isprime[n] <- !isprime[n] x <- x + 2 } y <- y + 6 } z <- z + 2 } y <- 3 while(y <= sqrt(limit)){ z <- 1 while(z <= 2){ x <- z while(x <= sqrt(limit) && (n <- 4*x*x+y*y) <= limit){ isprime[n] <- !isprime[n] x <- x + 3 } z <- z + 1 } y <- y + 6 } n <- 5 while(n <= sqrt(limit)){ if(isprime[n]){ k <- n*n while(k <= limit){ isprime[k] <- FALSE k <- k + n*n } } n <- n + 1 } isprime[2] <- isprime[3] <- TRUE if(return) (1:limit)[isprime] }
1000万までやろうとすると1分近くかかる。ダメっす。
比べてみよう
1000から1000万までの範囲でそれぞれの方法の計算速度を比べてみるとこんな感じ。単位は秒。
| range | sieve1 | sieve2 | sieve3 | sieve4 | atkin |
|---|---|---|---|---|---|
| 1e+03 | 0.001733333 | 0.001833333 | 0.001566667 | 0.003066667 | 0.003933333 |
| 1e+04 | 0.004666667 | 0.004133333 | 0.003866667 | 0.007566667 | 0.038533333 |
| 1e+05 | 0.022966667 | 0.016766667 | 0.013300000 | 0.022833333 | 0.456033333 |
| 1e+06 | 0.289600000 | 0.184900000 | 0.120100000 | 0.136800000 | 4.610400000 |
| 1e+07 | 2.974666667 | 1.869333333 | 1.281666667 | 1.275000000 | 44.410333333 |
グラフだとこんな感じ。
2の倍数、3の倍数、5の倍数...と飛ばす数を増やしていければ早いし探索範囲も広げられるんだろうけど自分の能力では効率よく実装できなくてなんともかんとも。
ちなみにグラフはこんな感じで書きました。
## グラフ化 r <- c(1e3, 1e4, 1e5, 1e6, 1e7) rep <- c(30, 30, 30, 10, 3) fname <- c("sieve1", "sieve2", "sieve3", "sieve4", "atkin") mean.elap.time <- function(f, index){ mean(sapply(1:rep[index], function(i){ system.time(eval(call(f, r[index])))[3] })) } elap.time.list <- function(f){ sapply(1:length(r), function(i){ mean.elap.time(f, i) }) } d <- data.frame( range = r, sieve1 = elap.time.list("sieve1"), sieve2 = elap.time.list("sieve2"), sieve3 = elap.time.list("sieve3"), sieve4 = elap.time.list("sieve4")) d2 <- data.frame(d, atkin = elap.time.list("atkin")) matplot(d[, 2:5], xaxt="n", xlab = "Upper limit", ylab = "Elapsed time (sec.)", pch = 1, type="l") axis(1, at=1:5, labels = r) legend("topleft", box.lty=0, legend=c("sieve1", "sieve2", "sieve3", "sieve4"), lty = 1:4, col = 1:4)
