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

10 ヶ月ぶりにシリーズ 6 回目の更新です。今の所これで打ち止め。

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

数学？

今回の高速化は第 3 回に続き 2 度目の、「計算理論」ではなく「計算機理論」から見たプログラム技法の紹介になります。

理論

各種パフォーマンス計測をしてみるとわかりますが、計算時間の大半を占めているのは SegmentedSieve() 内の以下の部分((3)より引用)になります。

while (j < sflags_.size()) {
  sflags_[j] &= kMask[ibit][k];
  j += i * C1[k] + C0[ibit][k];
  k = (k + 1) & 7;
}

ここで何をしているかというと

1. j < size の確認と分岐
2. sflags を更新
3. j の位置を更新
4. k を更新

です。ここで注意したいポイントは

• k は 0,1,2,3,4,5,6,7 の 8 値しか取り得ない
• j を進める長さは素数(i, ibit)を固定して考えると k 依存なので 8 パターンしかない
• j を更新する度に size との比較をしている

ということです。特に最初の k が 8 値しか取らないことはコンパイラ的には自明ではないためコンパイラの最適化で高速化されることはあまり期待できません。そこで人手で最適化を行います。

ループアンローリング (Loop unrolling)

今回紹介するループアンローリングはその他いろいろなプログラムでも応用できる手法です¹。唐突ですが、先のループにおいて、もし必ず k=0 からループスタートするとなれば以下のように書いても本質的に同じコードになるはずです。

k = 0;
while (j < sflags_.size()) {
  sflags_[j] &= kMask[ibit][0]; j += i * C1[0] + C0[ibit][0];
  if (j >= sflags_.size()) { k = 1; break; }  // (*)
  sflags_[j] &= kMask[ibit][1]; j += i * C1[1] + C0[ibit][1];
  if (j >= sflags_.size()) { k = 2; break; }  // (*)
  sflags_[j] &= kMask[ibit][2]; j += i * C1[2] + C0[ibit][2];
  if (j >= sflags_.size()) { k = 3; break; }  // (*)
  sflags_[j] &= kMask[ibit][3]; j += i * C1[3] + C0[ibit][3];
  if (j >= sflags_.size()) { k = 4; break; }  // (*)
  sflags_[j] &= kMask[ibit][4]; j += i * C1[4] + C0[ibit][4];
  if (j >= sflags_.size()) { k = 5; break; }  // (*)
  sflags_[j] &= kMask[ibit][5]; j += i * C1[5] + C0[ibit][5];
  if (j >= sflags_.size()) { k = 6; break; }  // (*)
  sflags_[j] &= kMask[ibit][6]; j += i * C1[6] + C0[ibit][6];
  if (j >= sflags_.size()) { k = 7; break; }  // (**)
  sflags_[j] &= kMask[ibit][7]; j += i * C1[7] + C0[ibit][7];
}

要するにループ内の k のパターンを全て列挙した形になります。この形でも少しコンパイラ最適化がかかりやすくなったかもしれませんが、もう 1 つ

• k=7 で sflags の更新を行うまで j < size となる

という仮定を置きましょう。つまり上記コード内 (**) の if 文必ず false になるという仮定です。すると ループ内で j < size の確認を行う演算と分岐 (上記コード内 (*) と (**)) が無駄になるので省略でき、何となく以下のようになりそうだということが分かります。

while (/* 何か条件 */) {
  sflags_[j] &= kMask[ibit][0]; j += i * C1[0] + C0[ibit][0];
  sflags_[j] &= kMask[ibit][1]; j += i * C1[1] + C0[ibit][1];
  sflags_[j] &= kMask[ibit][2]; j += i * C1[2] + C0[ibit][2];
  sflags_[j] &= kMask[ibit][3]; j += i * C1[3] + C0[ibit][3];
  sflags_[j] &= kMask[ibit][4]; j += i * C1[4] + C0[ibit][4];
  sflags_[j] &= kMask[ibit][5]; j += i * C1[5] + C0[ibit][5];
  sflags_[j] &= kMask[ibit][6]; j += i * C1[6] + C0[ibit][6];
  // assert (j < sflags_.size());
  sflags_[j] &= kMask[ibit][7]; j += i * C1[7] + C0[ibit][7];
}

見た目がすっきりしましたね。この形になると演算が減ったという点だけではなく、作業のパイプライン化による高速化も期待されますというかそちらの方が効果が高いです。

さて、次はその仮定を保証する「何か」を上手く入れていきましょう。

前処理

まずは前処理を行ってメインループに入る前に k=0 への調整を行います。これまでと同様の処理を回しつつ、k==0 になったらループを抜けるだけです。各 k について高々 1 回ずつしか回らないので、switch 文で書くこともできます。

// 前処理 (基本処理を k == 0 になるまで行う)
while (k && j < sflags_.size()) {
  sflags_[j] &= kMask[ibit][k];
  j += i * C1[k] + C0[ibit][k];
  k = (k + 1) & 7;
}

メインループ

計算コアの部分は while の条件をうまく調整することでコメントアウトしている asssert をクリアするようにしましょう。この assert の時点での j の値は、ループ先頭での j の値を $j_0$ とすると

\begin{eqnarray}
&&j_0 + \sum_{k=0}^{6} (i * C_1[k] + C_0[ibit][k])\\
&=&j_0 + 28i + S_{ibit}
\end{eqnarray}

となることがわかります。なお $S_{ibit}$ については

S_{ibit} = \sum_{k=0}^6 C_0[ibit][k] = m_{ibit} - 1

となるので、メインループの while 部分を以下のように変更することでループ内の分岐を消しても大丈夫ということができます。

while (j + i * 28 + kMod30[ibit] - 1 < sflags_.size()) {
  ...
}

後処理

あとは後処理ですが、必要なのは unroll したループの残り部分を処理することなので、実は元のルーチンがそのまま利用できます。

// 後処理 (元のままで実は問題ない)
while (j < sflags_.size()) {
  sflags_[j] &= kMask[ibit][k];
  j += i * C1[k] + C0[ibit][k];
  k = (k + 1) & 7;
}

ということで 前処理、メインルーチン、後処理を繋げたものを過去のメインルーチンと入れ替えることで高速化されるはずです。実際に変更を行ったコードがバージョン6 です。

計算時間

また毎度のごとく、計算時間を記録しておきます。今回の作業は前処理/後処理を追加したため $x\leqq 10^7$ では遅くなっていますが、$x$ が大きくなると高速になり $x=10^{10}$ で約 20% 高速化しています。

$x$	$\pi(x)$	Ver. 0	Ver. 1	Ver. 2	Ver. 3	Ver. 4	Ver. 5	Ver. 6
$10^7$	664579	0.034	0.011	0.004	0.004	0.004	0.012	0.015
$10^8$	5761455	0.437	0.148	0.060	0.056	0.052	0.061	0.036
$10^9$	50847534	-	3.157	1.652	0.576	0.593	0.460	0.348
$10^{10}$	455052511	-	-	19.806	6.347	6.351	5.013	3.766

区間篩の計算時間では大きな違いが見られません。これについては次回考察しましょう。

$x$	$y-x$	$\pi(y)-\pi(x)$	Ver. 4	Ver. 5	Ver. 6
$10^{12}$	$10^9$	36190991	0.782	0.656	0.553
$10^{13}$	$10^9$	33405006	0.942	0.760	0.683
$10^{14}$	$10^9$	31019409	1.197	1.129	0.963
$10^{15}$	$10^9$	28946421	1.948	1.825	1.817

1. ちなみに一般的なループアンローリングの手順はもう少し簡単です。 ↩