Search on the blog

2010年9月24日金曜日

エラトステネスの脅威

今日は、素数を列挙するアルゴリズムを考えます。
最初は基本的な実装を行い、徐々に高速化していきます。最終的には『エラトステネスの篩』を実装します。
また、それぞれのアルゴリズムのオーダーを見積り、実際に実行速度を計測します。

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



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]
getPrime1n^2157.9
getPrime2n^1.56.943
getPrime3n0.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)となります。

参考ページ:

0 件のコメント:

コメントを投稿