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

問題設定や以前の高速化については以下を参照してください。

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

省メモリ化

タイトルは高速化とありますが、今回の話題は省メモリ化です。
前回の高速化では一度に篩にかけるメモリ領域を絞りましたが、それはつまり篩ってない部分のメモリ情報は無くても問題無いということを意味しています。
なので、篩うためのメモリ領域を一度だけ確保し、そこを使いまわす形にしましょう。篩うことだけ考えると前回のメインルーチンである SegmentedSieve() は何も変えなくても問題ありません。考える必要があるのはそれを呼び出す元でメモリ領域の初期化の方法になります。

メモリ領域の初期化

では具体的にどのような注意が必要かというと、

• 1 に該当するビットを落とす
• 指定された上限 $x$ 以上に該当するビットを落とす

ということになります。(ここでの引用はバイトサイズの調整部分だけ)

void Sieve(const int64 x) {
  InitSflags(std::ceil(std::sqrt(x)));

  g_flags.resize(kSegmentSize, 0xff);
  for (int64 offset = 0, size = x / 30;  size > 0; offset += kSegmentSize, size -= kSegmentSize) {
    // サイズ調整
    if (size < kSegmentSize)
      g_flags.resize(size);

    if (offset)
      fill(g_flags.begin(), g_flags.end(), 0xff);
    else
      g_flags[0] = 0xfe;  // 1 のビットを削る
    SegmentedSieve(g_flags.data(), g_flags.size());
  }
}

省メモリ化のメリット

ところで、この省メモリ化を行うと微妙な計算や分岐が増えるので計算自体はむしろ少し（計測誤差の程度とはいえ）遅くなります。わざわざそんな変更を行ったメリットとしては、単純により広い領域の計算ができるようになることに加え、数値が大きな領域の部分的な篩い (例えば $[10^{17},10^{17}+10^{10})$ みたいな範囲の素数計算) ができるようになります。なので素数を求めていく作業を並列化、分散化することが可能になります。

範囲篩

ということで今回の後半は範囲篩 (Segmented Sieve) を行う際の変更点のお話です。以後、$[x, y)$ の範囲を篩にかけるという前提での変数設定です。

初期位置

素数 $p <\sqrt{y}$ で篩をかけるとき、通常では $p^2$ に該当するビットが篩の初期位置ですが、$p^2 < x$ となってしまっては計算ができません。一応、通常の篩を行うのと同じ手順で範囲内に篩位置が来るまですすめることも可能ですが、7 や 11 といった小さな素数に対して　$x = 10^{15}$ というような数が来る場合のことを考えるととても現実的ではありません。
なので、基本原理に立ち返りましょう。そもそも素数 $p=30k_1+m_1$ で篩っているとき、篩い落とされる合成数 $30k_c+m_c$ は

30k_c+m_c = (30k_1+m_1)(30k_2+m_2) = 30 (p k_2 + k_1 m_2) + m_1 m_2

という形で表されます。その中でも $x \leq 30k_c+m_c$ を満たす最小の $(k_2,m_2)$ を見つけたいのですが、直接計算で求めるのは無理そうです。なので

1. $m_2=1$ を仮定して $k_c \leq \lfloor x/30 \rfloor$ を満たす最大の $k_2$ を求める
2. 通常の篩と同じように篩位置を進める
3. $k_c \geq \lfloor x/30 \rfloor$ となったところから実際に篩を開始する

という手順で進めます。まずは $m_2=1$ を仮定すると初期位置は

30k_c+m_c = 30 (pk_2 + k_1) + m_1

となるので、$pk_2+k_1 \leq \lfloor x/30 \rfloor$ から

k_2 \leq \frac{1}{p} \left( \left\lfloor \dfrac{x}{30} \right\rfloor - k_1 \right)

で $k_2$ の初期値が求まります。コードに落とすと

const uint64 offset = x / 30;
for (uint64 i = 0; i < sflags_.size(); ++i) {
  for (uint8 primes = sflags[i]; primes; primes &= primes - 1) {
    uint8 lsb = primes & (-primes);
    int64 ibit = BitToIndex(lsb);
    const int64 m = kMod30[ibit];
    uint64 j = i * (30 * i + 2 * m) + (m * m) / 30;  // p^2 の場合の篩位置
    int64 jbit = ibit;
    if (j < offset) {
      int64 p = 30 * i + m;
      j = (offset - i) / p * p + i;  // i が数式の k_1、j が数式の k_c に該当する
      jbit = 0;
      while (j < offset) {
        j += i * C1[jbit] + C0[ibit][jbit];
        jbit = (jbit + 1) & 7;
      }
    }
    indecies.push_back(((j - offset) << 3) | jbit);
  }
}

ということになります。

計算時間

目安として、計算時間を記録しておきます。単位は秒。今回の更新は時間が短くなる期待はありませんが一応形式として。
今回のコードは github の Version 4 に該当します。

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

あとは区間篩の計算時間も。

$x$	$y-x$	$\pi(y)-\pi(x)$	Ver. 4
$10^{12}$	$10^9$	36190991	1.917