More than 5 years have passed since last update.

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

Last updated at 2019-07-11Posted at 2018-07-05

久々にシリーズ 5 回目の更新です。

初期フラグの設定

前回の変更で使うメモリ量が計算量的には $O(\sqrt{x})$ 程度まで、実際的にも MB 単位まで落ちてきたと思います。今回はその空いてできたスペースを使って行う高速化の紹介です。

理論

今回の高速化では 2 回目と同様に小さな素数での篩を省略します。
そもそもエラトステネスの篩を行う際、最初に全ビットを立てたフラグセットを用意しますが、この準備の際に全ビットを立てる必要はありません。合成数とわかっているものについては False になっていて構わないのです。
そのために「初期ビット列」のテンプレートを用意し、それをコピーして篩にかけることにします。テンプレートを作る際に使う素数 $p$ は面倒なので $p<30$ ということにしましょう。つまり

p_i = \{7, 11, 13, 17, 19, 23, 29\} = m_i

(の中のいくつか)ということになります。次はどの程度の長さでフラグの状態がループするのか考えてみましょう。元々 2, 3, 5 で篩った結果は ${\rm LCM}(2,3,5)=30$ を周期に $\phi(30)=8$ 個の候補があるのでそれを示すフラグをビット化していました。次に 7 で篩った状態を考えると、${\rm LCM}(30,7)=210$ を周期にループするはずです。今回はビットを詰める操作はしないので、ビット列的には 7 バイト¹を周期として $\phi(210)=48$ のビットが立っていることになります。同様に 11 でも篩うと 77 バイトで、13 でも篩うと 1001 バイトでループ、となっています。なのでその長さのビット列 (0xff で埋める) を作って、必要な素数で篩って初期ビット列を作ることにしましょう。

実装

実際には例えば 1001 バイトの 0xff 列を作って篩うのではなく、

2,3,5 で篩った結果として 1 バイトの 0xff がある
$p \in {7, 11, 13, 17, 19, 23}$ として
今の状態のビット列を $p$ 回繰り返す形に延長＆コピーする
そのビット列を $p$ で篩う

という手順にします。どうせ同じ列が続くのがわかっているのでコピーする方が楽です。あとは後の作業で楽になるので、本来の区間篩で使うサイズ $S$ だけの長さを先頭部分からコピーして後ろにくっつけておきましょう。
ちなみにこの初期ビット列を作るのに使う素数が 23 までなのは

7\cdot11\cdot13\cdot17\cdot19\cdot23=7436429

なので約 7 MB で済むのに対し、次の 29 も導入すると $7436429\cdot29=215656441$ と 200 MB になるからです。特にこのサイズで問題なければ使ってもらって構いません。

const int64 kInitialSize = 7 * 11 * 13 * 17 * 19 * 23;

void MakeInitialPattern() {
  int64 size = 1;
  g_initial.resize(size);
  g_initial[0] = 0xff;
  // {7, 11, 13, 17, 19, 23} に該当するビット位置をループ
  for (int64 ibit : {1, 2, 3, 4, 5, 6}) {
    const int64 p = kMod30[ibit];
    g_initial.resize(size * p);
    auto head = g_initial.begin();
    for (int64 i = 1; i < p; ++i)
      copy(head, head + size, head + size * i);
    size *= p;
    for (int64 j = 0, k = 0; j < size;) {
      g_initial[j] &= kMask[ibit][k];
      j += C0[ibit][k];
      k = (k + 1) & 7;
    }
  }
  // size == kInitialSize
  g_initial.resize(size + kSegmentSize);
  auto head = g_initial.begin();
  copy(head, head + kSegmentSize, head + size);
}

区間篩

これで区間篩本体の実装も少し変更する必要があります。具体的には

ビット列を 0xff で埋める形から、初期ビット列の該当部分をコピーする形に
篩う素数のリストから初期ビット列の形成に使ったものを除く

という変更になります。枠外チェックの部分を除いたコアな部分の実装だけ持ってくると

void Sieve(const int64 x) {
  InitSflags(std::ceil(std::sqrt(x)));
  sflags[0] = 0x80;  // 初期ビット列を作るのに使った素数を無くす

  g_flags.resize(kSegmentSize, 0xff);
  for (int64 offset = 0, size = x / 30;  size > 0; offset += kSegmentSize, size -= kSegmentSize) {
    // 初期ビット列のコピー
    auto head = g_initial.begin() + offset % kInitialSize;
    copy(head, head + kSegmentSize, g_flags.begin());
    if (offset == 0)
      g_flags[0] &= 0xfe;

    SegmentedSieve(g_flags.data(), g_flags.size());
  }
}

という形になります。SegmentedSieve() 内で使う sflags の変更を関数冒頭で行っているので見落とさないようにしましょう。

計算時間

また毎度のごとく、計算時間を記録しておきます。単位は秒。今回は時間計測を行っている VM とコンパイラが新しくなったおかげで以前の記事と比較すると軒並み倍速くらいになっていますので、この表の中だけで比較してください。
今回のコードは github の Version 5 に該当します。

$x$	$\pi(x)$	Ver. 0	Ver. 1	Ver. 2	Ver. 3	Ver. 4	Ver. 5
$10^7$	664579	0.023	0.010	0.004	0.004	0.004	0.003
$10^8$	5761455	0.375	0.121	0.058	0.047	0.050	0.036
$10^9$	50847534	-	2.831	1.425	0.540	0.572	0.425
$10^{10}$	455052511	-	-	18.685	5.846	6.190	4.822

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

$x$	$y-x$	$\pi(y)-\pi(x)$	Ver. 4	Ver. 5
$10^{12}$	$10^9$	36190991	0.756	0.614
$10^{13}$	$10^9$	33405006	0.880	0.717
$10^{14}$	$10^9$	31019409	1.111	0.987
$10^{15}$	$10^9$	28946421	1.741	1.560

1 バイトが 8 ビットでないことをは想定してません ↩

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up