エラトステネスの脅威

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

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

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

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

view plainprint?

1. void getPrime1(int n) {  
2.     VI prime;  
3.   
4.     FOR (i, 2, n) {  
5.         bool isPrime = true;  
6.   
7.         FOR (j, 2, i) {  
8.             if (i % j == 0) {  
9.                 isPrime = false;  
10.                 break;  
11.             }  
12.         }  
13.         if (isPrime)  
14.             prime.PB(i);  
15.     }  
16.   
17.     printf("%d\n", prime.size());  
18. }  

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)までの割り算をすればよいことに気付く。これをふまえて、改良するとこうなる。

view plainprint?

1. void getPrime2(int n) {  
2.     VI prime;  
3.   
4.     FOR (i, 2, n) {  
5.         bool isPrime = true;  
6.   
7.         FOR (j, 2, (int)sqrt(i) + 1) {  
8.             if (i % j == 0) {  
9.                 isPrime = false;  
10.                 break;  
11.             }  
12.         }  
13.         if (isPrime)  
14.             prime.PB(i);  
15.     }  
16.   
17.     printf("%d\n", prime.size());  
18. }  

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());
}

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

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

view plainprint?

1. void getPrime3(int n) {  
2.     bool *isPrime = new bool[n];  
3.     VI prime;  
4.   
5.     memset(isPrime, 1, sizeof(bool) * n);  
6.     for (int i = 2; i < sqrt(n) + 1; i++) {  
7.         if (!isPrime[i])  
8.             continue;  
9.   
10.         for (int j = i + i; j < n; j += i)  
11.             isPrime[j] = 0;  
12.     }  
13.     FOR(i, 2, n)  
14.         if (isPrime[i])  
15.             prime.PB(i);  
16.   
17.     delete [] isPrime;  
18.     printf("%d\n", prime.size());  
19. }  

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のときの実行速度を計測してみましょう。

view plainprint?

1. #include <sys/time.h>  
2.   
3. inline double gettimeofday_sec() {  
4.     struct timeval t;  
5.     gettimeofday(&t, NULL);  
6.     return (double) t.tv_sec + (double) t.tv_usec * 1e-6;  
7. }  
8.   
9. int main() {  
10.     double s = gettimeofday_sec();  
11.     getPrime1(1000000);  
12.     printf("%lf\n", gettimeofday_sec() - s);  
13.   
14.     s = gettimeofday_sec();  
15.     getPrime2(1000000);  
16.     printf("%lf\n", gettimeofday_sec() - s);  
17.   
18.     s = gettimeofday_sec();  
19.     getPrime3(1000000);  
20.     printf("%lf\n", gettimeofday_sec() - s);  
21.   
22.     return 0;  
23. }  

#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