エラトステネスの篩

正の整数N未満の素数を求める

プログラム

C

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
#include <math.h>
#include <stdio.h>
#include <time.h>

#define N 100000000
unsigned flags[N / 8];
#define INT_SIZE 32
#define DOUBLE_INT_SIZE 64
#define GETBIT(n) (flags[n / DOUBLE_INT_SIZE] & 1 << (n / 2 % INT_SIZE))
#define SETBIT(n) (flags[n / DOUBLE_INT_SIZE] |= 1 << (n / 2 % INT_SIZE))
#define UNSETBIT(n) (flags[n / DOUBLE_INT_SIZE] ^= 1 << (n / 2 % INT_SIZE))

void eratosthenes(int n) {
    for (int i = 3, m = floor(sqrt(n)); i < m; i += 2) {
        if (GETBIT(i)) continue;
        for (int j = i * i; j < n; j += 2 * i) SETBIT(j);
    }
}

void printflags(int max) {
    for (int i = 1, c = 0; i < max; i += 1) {
        if (i % 2 ? !GETBIT(i) : i == 2) {
            printf("%d, ", i), c++;
            if (!(c % 20))
                putc('\n', stdout);
        }
    }
}

int main() {
    struct timespec f, t;
    timespec_get(&f, TIME_UTC);

    eratosthenes(N);

    timespec_get(&t, TIME_UTC);
    long nsec;
    long long tsec;
    if (f.tv_nsec < t.tv_nsec) {
        nsec = t.tv_nsec - f.tv_nsec;
        tsec = t.tv_sec - f.tv_sec;
    } else {
        nsec = 1000000000 + t.tv_nsec - f.tv_nsec;
        tsec = t.tv_sec - f.tv_sec - 1;
    }
    printf("%lld.%09ld\n", tsec, nsec);
    // printflags(N);
}
1
g++ eratosthenes.c -o eratosthenes.out -lm -O3

Ruby

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
def eratosthenes(n)
  flags = Array.new(n, 0)
  start = Time.now()
  3.step(Math::sqrt(n).floor(), 2) do |i|
    if flags[i] == 1
      next
    end
    (i * i).step(n, i * 2) do |j|
      flags[j] = 1
    end
  end
  p Time.now() - start
  1.step(n, 1) do |i|
    if (i % 2 == 1) ? (flags[i] == 0) : (i == 2)
      printf("%d, ", i)
    end
  end
end

eratosthenes(100)
# 1.08e-05s
# 1, 2, 3, 5, 7, 11, 13, 17, 19, 23,
# 29, 31, 37, 41, 43, 47, 53, 59, 61,
# 67, 71, 73, 79, 83, 89, 97, 
eratosthenes(100000000)
# 3.804s (表示は省略)

結果

秒/N 1000万 1億
C 0.007 0.088
Ruby 0.346 3.804

(Cの最適化は-O3)

考察

1
2
3
4
5
6
void eratosthenes(int n) {
    for (int i = 3, m = floor(sqrt(n)); i < m; i += 2) {
        if (GETBIT(i)) continue;
        for (int j = i * i; j < n; j += 2 * i) SETBIT(j);
    }
}

アルゴリズム

1. (全体)2の倍数で素数なのは2だけだから、forを2の倍数で回す。2は表示するときに例外的に処理する

2の倍数を考えないことで計算量減る。また、1と2を足して回るときに、機械語数は変わらないので性能の劣化はない。

2. (外側for)3から平方根未満の数に対して合成数をマークしていく

N未満の数のそれぞれの合成数の素因数の少なくとも１つはNの平方根未満であるからである。3から回すのは、1.で、2の倍数は考えないようにしているからである。

3. (内側for)i**2からループを始める

i*2の倍数でこれより小さいものはすでに除外済みである。i**2未満のそれぞれの合成数の素因数の少なくとも１つはi未満であるからである。

4. (外側for内のif)倍数で消していこうとしている数が合成数ならば、それの倍数は既に除外済みであるので無視する

iが今回のiの素因数のときに既に除外済み

5. forの範囲において、無駄な計算はしない

1
2
3
for (int i = 3; i < (int)floor(sqrt(n)); i += 2)
// 上のようにせずに下のようにする
for (int i = 3, m = floor(sqrt(n)); i < m; i += 2)

floor(sqrt(n))の計算結果はループ中に変わらない。

最適化

6. ビット単位で素数か合成数かを記録する

bit単位ではなく、int単位(4byte)、char単位でも計測してみた。
それぞれN=1億の場合である。bit単位で記録すればメモリの無駄は0であるが、除算・剰余演算が必要で計算量は増大する。対して、char、int単位で記録すればメモリの無駄はそれぞれ7bit*N、31bit*Nであるが、割り算は必要ない。

計測の結果は次のようになった。

記録単位 時間(秒)
bit 0.088
char 0.348
int 0.549

結果はbit単位が最も高速だった。このことから、bit単位のほうがCPUのキャッシュヒット率が高いため高速になっているものと思われる。CPUは非常に高速なのに対してメモリアクセスは非常に遅いので、多少CPUの計算量が増加しても、なるべくキャッシュにアクセスできたほうがよい事がわかる。

7. マクロを使い、関数呼び出しを少なくする

bitを立てるときに関数を呼び出すのと、マクロを使うのとでは速度に差が出た。

1
2
#define GETBIT(n) (flags[n / DOUBLE_INT_SIZE] & 1 << (n / 2 % INT_SIZE))
#define SETBIT(n) (flags[n / DOUBLE_INT_SIZE] |= 1 << (n / 2 % INT_SIZE))
1
2
3
4
5
6
int GETBIT(int n){
    return flags[n / DOUBLE_INT_SIZE] & 1 << (n / 2 % INT_SIZE);
}
void SETBIT(int n){
    flags[n / DOUBLE_INT_SIZE] |= 1 << (n / 2 % INT_SIZE);
}
時間(秒)
関数呼び出し(-O0) 0.313
マクロ使用 0.088

マクロを使用し、関数呼び出しを少なくすると、高速になった。このような小さな関数では機械語レベルで関数呼び出しコストが無視できなくなるからである。