エラトステネスの篩の高速化 (3)

エラトステネスの篩の高速化シリーズ3本目です。問題設定は 1 本目を参考に。
この辺りの高速化からは Qiita の他の記事に載ってない気がします。[要出典]

• エラトステネスの篩の高速化 (1)
• エラトステネスの篩の高速化 (2)
• エラトステネスの篩の高速化 (3) ← 今ココ
• エラトステネスの篩の高速化 (4)
• エラトステネスの篩の高速化 (5)
• エラトステネスの篩の高速化 (6)

使うメモリ容量

2,3,5 の倍数を除いた実装にするとフラグ管理に使うメモリ容量はだいたい $x/30$ バイトとなるので、例えば $x=10^9$ で 30MB、$x=10^{10}$ で 300MB のメモリを使うことになります。また同時に、それだけのメモリ領域全域に渡って読み書きを行うのでキャッシュミスがよく起きますというか大体キャッシュミスになります。
これを高速化するために、篩が一括で処理する範囲を小さくしてキャッシュメモリに乗るようにするのが今回のお話です。

区間篩 (Segmented Sieve)

ということで、適当なバイト数 $S$ ごとにメモリを区切って篩います。篩う区間ごとに素数をループさせますが、次の区間に移ると篩う位置を計算で求めることになります。これについては毎回計算すると面倒なのでメモしておくことにしましょう。

篩い位置のメモ

前回の処理に出てくる $(i_2,j_2)$ を素数と対応させて記録することにします。使う数の範囲としては $0 \leqq i_2 < x/30$、$0\leqq j_2 < 8$ なので 2 つまとめて

  int64 k2 = (i2 << 3) + j2;

として保持することにします。ちなみに篩に使う $\sqrt{x}$ までの素数については計算が重複するかもしれませんが、別の配列に管理しておきます。これによりメモリ量がフラグ管理で $\sqrt{x}/30$ バイト、上記の篩い位置の保存に $8\pi(\sqrt{x}) \simeq 16\sqrt{x}/\ln x$ バイト程度増えますが、篩本体の対象に比べれば小さいので気にしないことにしましょう。

基本的には入力になる $x$ に対して $\sqrt{x}$ までの素数とそれに対応する篩い位置がリストアップされればいいので、作業量を最小化するには $\sqrt[4]{x}$ までの素数で篩いつつ位置計算をし、それ以上の数については位置計算をするだけで良いのですが、処理が面倒なのでここでは $\sqrt{x}$ (= 入力になる x) まで篩う形で書いてます。

std::vector<uint8> sflags;  // sqrt(x) までの素数フラグ。
std::vector<uint64> indecies;  // i 番目の素数での篩い位置

// x 以下の素数リストを作って篩位置を初期化する
void InitSflags(const int64 x) {
  const int64 size = x / 30;
  for (uint64 i = 0; i < size; ++i) {  // 本当は sqrt(x)/30 までで篩えば十分です。
    for (uint8 flags = sflags[i]; flags; flags &= flags - 1) {
      uint8 lsb = flags & (-flags);
      int ibit = BitToIndex(lsb);
      const int32 m = kMod30[ibit];
      int32 pm = 30 * i + 2 * m;
      uint64 j = i * pm + (m * m) / 30;
      uint64 k = ibit;
      indecies.push_back(((j + kSegmentSize) << 3) | k);
      while (j < sflags_.size()) {
        sflags_[j] &= kMask[ibit][k];
        j += i * C1[k] + C0[ibit][k];
        k = (k + 1) & 7;
      }
    }
  }
}

篩い本体

本来の篩を行う本体は先に述べたとおり、適当に決めたサイズごとにエリアを区切って篩います。区間を変えたときにその辺りの端数処理さえ気をつければ難しくないと思います。

const uint64 kSegmentSize = 2 * 1024 * 1024;  // 例えば 2MB ちょうど
std::vector<uint8> g_flags;  // 計算対象になるフラグ全部

void Sieve(const int64 x) {
  int64 size = x / 30;
  g_flags.resize(size, 0xff);
  g_flags[0] = 0xfe;
  InitSflags(std::ceil(std::sqrt(x)));
  for (uint8* flags = g_flags.data(); size > 0; flags += kSegmentSize, size -= kSegmentSize) {
    SegmentedSieve(flags, size);
  }
}

void SegmentedSieve(uint8* flags, const int64 size) {
  auto p_index = indecies.begin();
  for (uint64 i = 0; i < sflags.size(); ++i) {
    for (uint8 primes = sflags[i]; primes; primes &= primes - 1) {
      uint8 lsb = primes & (-primes);
      int ibit = BitToIndex(lsb);
      uint64 index = *p_index;
      uint64 j = (index >> 3) - kSegmentSize;
      uint64 k = index & 7;
      while (j < size) {
        flags[j] &= kMask[ibit][k];
        j += i * C1[k] + C0[ibit][k];
        k = (k + 1) & 7;
      }
      *p_index = (j << 3) | k;
      ++p_index;
    }
  }
}

計算時間

目安として、計算時間を記録しておきます。単位は秒。今回は篩う素数リストを作るために計算量が少し増えた一方、$x\geqq 10^8$ では本篩部分がインキャッシュになって高速化したので全体としては高速になってます。
今回のコードは github の Version 3 に該当します。

$x$	$\pi(x)$	Ver. 0	Ver. 1	Ver. 2	Ver. 3
$10^7$	664579	0.047	0.017	0.007	0.006
$10^8$	5761455	1.039	0.352	0.189	0.076
$10^9$	50847534	-	6.180	3.377	0.898

次回 は今回の応用として、利用するメモリ量を減らす方法を紹介します。