パソコンなど

パソコンの話をします。

エラトステネスの篩の高速化

　少し前に某会社の採用試験問題(?)の解説みたいな感じでエラトステネスの篩が話題になっていたので書こうと思う。C++を使うよ。

エラトステネスの篩とは

  • ある数以下の素数を比較的高速に求めるアルゴリズム。数の集合から素数の倍数を取り除くことで最後に素数が残る。
  • 以下はlimit以下の素数を格納したstd::vector<int>を返す関数の実装。
std::vector<int> simple_eratosthenes(const long long limit) {
    std::vector<bool> data(limit+1, false);
    data[0] = true;
    data[1] = true;
    for (long long s = 2; s*s <= data.size();) {
        for (long long i = s*s; i < data.size(); i+=s) {
            data[i] = true;
        }
        s++;
        while (data[s] == true) {
            s++;
        }
    }
    std::vector<int> result;
    for (long long i = 0; i < data.size(); i++) {
        if (data[i] == false) {
            result.push_back(i);
        }
    }
    return result;
}

やったこと

  • 最初から偶数を取り除いておくのと同様に、2, 3, 5, 7...など任意の数までの倍数を取り除けるようにした。
  • 2の場合だと、 2x+1のみが残る。
  • 3を含めると、 6x+1, 6x+5のみが残る。
  • 5を含めると、 30x+1, 30x+7, 30x+11, 30x+13, 30x+17, 30x+19, 30x+23, 30x+29 のみが残る。
  • これらのそれぞれを別々の実行スレッドに投げて並列化する。
  • 2の倍数を取り除くと0.5倍、3の倍数まで取り除くと0.33倍、5の倍数まで取り除くと0.267倍というようにメモリの消費量も減る。
  • それぞれに素数は同じ感じで分布しているらしいので並列化のムラはほとんどない。

実装関係の話

  • 上記のsimple_eratosthenesを使って探索範囲の平方根までの素数を求める。
    • simple_eratosthenesはメインの探索を行いながら、動的に篩うための素数を探しているが、マルチスレッド化するためにはこの方法だと別のスレッドからデータを読まなければならない。これはだるい。シングルで事前に計算した方がいい。
  • 各スレッドで、各素数を篩うときに開始位置を算出する必要がある。(例えば、 30x+13=41yを満たすxを求める。(答えは 943=41\times 23=30 \times 31+13\therefore x=31 ))
    • これは二元一次不定方程式なので、拡張互除法などで比較的高速に解ける。
    • 剰余が1となる場合だけ求めておけばあとはそれを何倍かすれば他はもとまるので、スレッドを分割する前に1についてのみを計算しておき、あとはスレッドに任せる。
    • 例えば、 30x_1+1=41y_1を満たす最小の自然数 x_1 は、 451 = 30\times15+1 = 41\times11\therefore x_1=15となる。
    • 30x_2+13=41y_2を満たす自然数 x_2 は、 x_2\equiv x_1\times13\bmod41である。 x_2\equiv195\bmod41なので x_2=31となる。
  • シングルで求める比較的計算量が大きい部分はこの二つだが計算量は O(\sqrt{n \log n})程度なので、メモリの確保やメインの処理に比べれば雀の涙。(メモリの確保の方が圧倒的に時間がかかる。)

時間計測

  • 私のノートパソコン(MacBook Pro (13-inch, 2018, Four Thunderbolt 3 Ports))で行う。CPUはIntel(R) Core(TM) i5-8259U CPU @ 2.30GHzで、メモリは8GB 2133 MHz LPDDR3
  • コンパイラはgcc version 10.2.0 (Homebrew GCC 10.2.0_4)を使い、最適化オプションを盛る。
  • 前述の篩うための素数を探す計算と、開始位置を探す計算は無視する。(1012までやっても14msしかかからない。)
  • メモリの確保にかかる時間は無視する。(今回やりたい評価には無関係なので)
  • 結果の素数をカウントし、その個数を以て正しさを評価する。
  • 素数の個数については、みんなの知識 ちょっと便利帳に掲載されているものを使う。

結果

  • 何個の素数を使って分割するかをbase_sieve_sizeで表す。
  • thread数、圧縮率以外の単位はミリ秒
  • すべて一発取りなので結果はあくまで目安
base_sieve_size 3 4 5 6 7
thread数 8 48 480 5760 92160
圧縮率 3.75 4.375 4.8125 5.21354 5.53939
107 1 1 12 309 44406
108 10 6 13 307 52329
109 684 109 57 325 60774
1010 15246 7033 907 696 62785
1011 208788 177030 68373 8994 58088
  • 結果は全て正しかった。

考察など

  • 問題サイズを大きくした場合、大量のスレッドに分割した方が速くなっている。
    • 問題を細分化したことにより、一つ一つのスレッドがキャシュに乗るようになり、動作が高速になったものと思われる。
    • base_sieve_size = 7の時は基本的に遅いが、問題サイズによる変化が最も小さいので、1012以上の大きさならばより高速に動作する可能性がある。
  • 1011をbase_sieve_size = 6で計算するために、すでに2.23GBのメモリ領域を使っている。メモリが足りないので1012までを一気に求めることはできない。
    • 計算結果を順次捨てながらメモリ領域を次のスレッドに引き継ぐようにすればメモリ不足は解消できるが、それは別の機会にやると思う。

ソースコード

github.com