Project Eulerの素数を扱う問題で必要になったので、エラトステネスの篩を実装してみた。
実装した関数sieve
は、上限max
までの数の中で素数であるものを探し、対応するIsPrime
の要素をtrue
にする。
#うろ覚え
とりあえず何も見ずに実装してみたら、違うものになってしまった。
std::vector<bool> IsPrime;
void sieve(size_t max){
if(max+1 > IsPrime.size()){ // resizeで要素数が減らないように
IsPrime.resize(max+1,true); // IsPrimeに必要な要素数を確保
}
IsPrime[0] = false; // 0は素数ではない
IsPrime[1] = false; // 1は素数ではない
for(size_t i=0; i<=max; ++i) // 0からmaxまで全て調べる
if(IsPrime[i]) // iが素数ならば
for(size_t j=i+1; j<=max; ++j) //i+1からmaxまで全て調べる
if(j%i==0) // iの倍数は
IsPrime[j] = false; // 素数ではない
}
最初はこれが正しいものと思い込んで使っていたが、あまりにも遅いので怪しく思い始めた。
##なぜ遅い?
・外側のループが無駄に回っている
合成数を除くためならi<=sqrt(max)
(i*i<=max
)までで十分
(合成数$n$は必ず$\sqrt{n}$以下に約数を持つ)
・内側のループが非効率
わざわざ全ての数について、i
の倍数かを調べている
最初からi
の倍数だけを回れば済む話
(思いがけずも)やっていることは素朴な試し割りと同じなので、遅いのも当然。
#ググった
ググってちゃんとした実装を書いた。
std::vector<bool> IsPrime;
void sieve(size_t max){
if(max+1 > IsPrime.size()){ // resizeで要素数が減らないように
IsPrime.resize(max+1,true); // IsPrimeに必要な要素数を確保
}
IsPrime[0] = false; // 0は素数ではない
IsPrime[1] = false; // 1は素数ではない
for(size_t i=2; i*i<=max; ++i) // 0からsqrt(max)まで調べる
if(IsPrime[i]) // iが素数ならば
for(size_t j=2; i*j<=max; ++j) // (max以下の)iの倍数は
IsPrime[i*j] = false; // 素数ではない
}