エラトステネスの篩 in C++

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;      // 素数ではない
}