エラトステネスの脅威

今日は、素数を列挙するアルゴリズムを考えます。

最初は基本的な実装を行い、徐々に高速化していきます。最終的には『エラトステネスの篩』を実装します。

また、それぞれのアルゴリズムのオーダーを見積り、実際に実行速度を計測します。

まずは、基本的な実装から。これ。

void getPrime1(int n) {
    VI prime;

    FOR (i, 2, n) {
        bool isPrime = true;

        FOR (j, 2, i) {
            if (i % j == 0) {
                isPrime = false;
                break;
            }
        }
        if (isPrime)
            prime.PB(i);
    }

    printf("%d\n", prime.size());
}

まー、だれでも思いつくような方法。

少し考えてみると、整数iが素数か否かを判定するのにi-1までの割り算をしなくても、sqrt(i)までの割り算をすればよいことに気付く。これをふまえて、改良するとこうなる。

void getPrime2(int n) {
    VI prime;

    FOR (i, 2, n) {
        bool isPrime = true;

        FOR (j, 2, (int)sqrt(i) + 1) {
            if (i % j == 0) {
                isPrime = false;
                break;
            }
        }
        if (isPrime)
            prime.PB(i);
    }

    printf("%d\n", prime.size());
}

ちょっと変わっただけですが、かなり速くなりそうです。

次に、「エラトステネスの篩」を用いた実装。（アルゴリズムの考え方については、参考ページ参照のこと。）ソースを見てもらえば何をしているか分かると思います。

void getPrime3(int n) {
    bool *isPrime = new bool[n];
    VI prime;

    memset(isPrime, 1, sizeof(bool) * n);
    for (int i = 2; i < sqrt(n) + 1; i++) {
        if (!isPrime[i])
            continue;

        for (int j = i + i; j < n; j += i)
            isPrime[j] = 0;
    }
    FOR(i, 2, n)
        if (isPrime[i])
            prime.PB(i);

    delete [] isPrime;
    printf("%d\n", prime.size());
}

ここまで、3つの方法を記述しました。それぞれの計算量を見積ってみましょう。そして実際にn=1,000,000のときの実行速度を計測してみましょう。

#include <sys/time.h>

inline double gettimeofday_sec() {
    struct timeval t;
    gettimeofday(&t, NULL);
    return (double) t.tv_sec + (double) t.tv_usec * 1e-6;
}

int main() {
    double s = gettimeofday_sec();
    getPrime1(1000000);
    printf("%lf\n", gettimeofday_sec() - s);

    s = gettimeofday_sec();
    getPrime2(1000000);
    printf("%lf\n", gettimeofday_sec() - s);

    s = gettimeofday_sec();
    getPrime3(1000000);
    printf("%lf\n", gettimeofday_sec() - s);

    return 0;
}

結果は以下のとおりになりました。

アルゴリズム	オーダー	実行時間[s]
getPrime1	n^2	157.9
getPrime2	n^1.5	6.943
getPrime3	n	0.016

やっぱりエラトステネス速いですね。。getPrime1, getPrime2のオーダーは見てのとおりですが、getPrime3（エラトステネス）のオーダーは少し分かり難いので少し補足します。

まず、一番外側のループは見てのとおりO(n^0.5)です。次に内側のループですが、n/i回まわります。

iによって値が違うので平均値を考えます。i は 2 から n^0.5 までなので、内側のループの回数は、それぞれのiに対して、n/2, n/3, n/4, ...... n/n^0.5となります。この平均は、およそn/{(2+n^0.5)/2}なので内側のループの回数はO(n^0.5)です。よって2重ループ部分のオーダーは、O(n^0.5) * O (n^0.5) = O(n)です。最後にisPrimeを線形探索する部分もO(n)なので、最終的にエラトステネスのオーダーは、O(n)となります。

参考ページ:

エラトステネスの篩について　http://ja.wikipedia.org/wiki/%E3%82%A8%E3%83%A9%E3%83%88%E3%82%B9%E3%83%86%E3%83%8D%E3%82%B9%E3%81%AE%E7%AF%A9