えびちゃんの日記

えびちゃん(競プロ)の日記です。

全方位木 DP を指して rerooting と呼ぶ風潮に対するお気持ち表明

rerooting DP、rerooting をする DP、DP with rerooting とかなら別によさそう。

全方位木 DP 自体を指して「rerooting」と言われると気になります。 rerooting は木の根を変える操作(全方位木の中で考えられる)を指していて、全方位木 DP 自体のことを指しているわけではないと思っています。

例を出すよ

DP 以外の文脈で rerooting という言葉が使われているのを挙げてみます。 Euler tour を管理するデータ構造で、根とする頂点を変える操作が rerooting と呼ばれています。

http://web.stanford.edu/class/archive/cs/cs166/cs166.1166/lectures/17/Slides17.pdf

上の資料では、森に関するクエリ処理をするデータ構造を考えています。 対象のクエリは「この辺を追加してね」「この辺を削除してね」で、ループはできないことが保証されています。

森を Euler tour の集合として表現すると、辺の追加や削除に際して、Euler tour の根の頂点を変える操作ができると楽になるようです。 そこで、その操作名が rerooting と呼ばれています。

例おわり。

DP の名前ではなくて、単に、より一般的な操作の名前として捉える方が自然な単語に思えます。

そういえば

splay tree ではアクセスしたノードを根に持っていく操作を splay と呼んでいましたね*1

そこで思ったんですが、splay tree で splay をする際に DP の更新をするようにすれば、順に頂点をなめることで全方位木 DP で得たいものを得られたりしませんかね? splay は amortized \(O(\log(n))\) time ですが、順になめたときは amortized \(O(\log(\log(n)))\) time という結果もあった気がします(詳細見失いました、嘘かも)。

何にせよ普通に BFS/DFS した方が速そうですが...

おわり

もしかしたらえびちゃんの勘違いで、全方位木 DP 自体を指して rerooting と呼んでいる人はいないかもしれません?

お気持ち表明シリーズが好きな人におすすめ(今回の記事より内容ずっしりめです): rsk0315.hatenablog.com

*1:link/cut tree では evert という呼ばれ方もあったような...?

インタラクティブ問題のデバッグに関してなど

インタラクティブ問題は手元でのデバッグが面倒がちなので、いやだと思っている人が多そうです。 coproc というシェルの機能を使えばそこそこ簡単にできないかな?と思ったので、書いてみます。

まえおき

簡単とは言っても、以下のものは自分で書く or 生成する必要があります*1

ジャッジについては、次のような仕様を満たしていてほしいです。

  • 入力データのファイル名をコマンド引数で受け取る
    • あるいは、ジャッジがランダムで生成してもよいです、ここはよしなに。
    • 標準入力はクエリとのやりとり用にしたいので使えません。
  • クエリを標準入力から受け取る
  • 返答を標準出力に書き出す
  • デバッグ用の出力を標準エラー出力に出してもよい
  • AC ならステータス 0 で終了、そうでなければ 1 などで終了*2
    • 要するに、return ac? 0: 1; みたいなことをしてほしい。

やります

さて、解答のプログラム ./sol とジャッジのプログラム ./judge と入力ファイル infile を用意できたとします。 このとき、次のようにできます。シェルによって違うので、各々説明します。

Bash

実行例:

bash-5.1$ coproc ./sol
[1] 27352
bash-5.1$ ./judge infile <&${COPROC[0]} >&${COPROC[1]}
AC: ok
[1]+  Done                    coproc COPROC ./sol

実行例なので出力も混ざっていますが、実際に書く必要があるのは次の部分だけです。

coproc ./sol
./judge infile <&${COPROC[0]} >&${COPROC[1]}

coproc をつけてコマンドを実行すると、標準入出力以外と入出力をやりとりするものとして開始されます(パイプと呼ばれる機構とやりとりします)。 入力先のパイプは file descriptor*3 の形で教えられ、その値は ${COPROC[1]} に入っています。出力先も同様で、${COPROC[0]} です。

file descriptor が fd のものにリダイレクトしたいときは、<&fd>&fd の形を使えばよいです。 標準エラー出力を標準出力にまとめるときに 2>&1 というのを見たことがある人もいると思いますが、標準出力の file descriptor が 1 なので、こういう感じの書き方になっているわけですね。 標準入力は 0標準エラー出力2 です。

coproc をつけて実行したコマンドにとっての入力先は ${COPROC[1]} に入っていますが、それに対して書き込む(入力を与える)コマンドにとっては ${COPROC[1]} は出力先なので、標準出力の番号に対応する 1 に値が入っていると考えるとよさそうです。若干ややこしいですね。

AC: ok標準エラー出力に対する書き込みです。実際には、次のような内容があると便利かもしれません。

  • どういうクエリを受け取ったか
  • どういう返答をしたか
  • QLE や WA であればその詳細

Zsh

実行例:

% coproc ./sol
[1] 27741

% ./judge infile <&p >&p
[1]  + done       ./sol
AC: ok

書く必要がある部分:

coproc ./sol
./judge infile <&p >&p

こっちの方が短いですね。Zsh では、coproc で作ったコマンドに対してリダイレクトしたいときは <&p>&p が使えるとのことなので、それを使えばよいです。 他は Bash と同様です。

そのほか

IDE とかを使っていてシェルと仲よくない人もいるらしいです? よしなにがんばってほしいです...

あと fish-shell の人もよしなにがんばってほしいです。

そうなんだけど

github.com

oj っていうツールがあるんですが、次のようなことができます。

% oj t/r -c ./sol './judge infile'

おまけ

macOS にデフォルトで入っている Bash のバージョンは化石であることで有名です。

bash-3.2$ coproc
bash: coproc: command not found

brew install bash とかをするとよいです。

おわり

結局、ジャッジを書くのが面倒という問題からは逃れられないので困ります(これは構築などのスペシャルジャッジが必要な問題についてもそう)。 ジャッジを書くパートさえできたら、解答とジャッジのやりとりパートはわりと簡単じゃんというお話でした。

ジャッジを書くパートは面倒で、特に adaptive/adversarial な問題では、厳しく書くのは大変です。 これは、クエリのやりとりに矛盾しない入力データであれば、ジャッジ側が "後出し" できるといった問題設定のことです。

それはそれとして、interactive や reactive などの表記揺れ(表記というか表現というか)があって、どれを使うべきなんだろう〜と思っています。

*1:やっぱり面倒じゃん、解散。←?

*2:ステータスを -1 などの負数にする人がよくいますが、これは予約された値であり、よくなかった気がします。文献は見失いました。

*3:ここでは、単にファイルやパイプを特定するための表し方の一つ、くらいに認識しておけばいい気がします。

誤差許容ジャッジちゃんと遊んでみる (part 0)

part 1 以降があるかはわかりません。

導入?

競プロで、「絶対誤差または相対誤差が \(10^{-6}\) 以下ならば正解と判定される」といった文言はちょくちょく見かけます。

思い出すシリーズ \(2e-3\) atcoder.jp

その昔、これに対して NaN と出力するだけで AC になるなどの脆弱性もありました。

atcoder.jp

最近の問題のジャッジではそうしたことは起こらないはずです。

atcoder.jp

自前で書くとそういう脆弱性を入れがちなので、testlib とかのライブラリなどを適切に使うとよいと思います。 普段 contestant として「お行儀のよい」入力に慣れがちな競プロ er は、そうした脆弱性に疎いんじゃないかな?という気がしています(偏見)。

rian.hatenablog.jp

本題?

contestant 的には、たとえば答えが 1.5 のときに 1.50 とか 1.50000000000000000 とか 1.4999999999999999 とかを出力しても WA になってほしくないですよね。 あるいは、指数表記を許すと言われているジャッジであれば、15e-1 とか 150000000000000000000e-20 とか 0.000000000149999999999999e+10 とかを出力しても許されたいわけです。

ということで、ジャッジちゃんと遊びました。

1.5000...00 の形式 atcoder.jp

15000...00e-xx の形式 atcoder.jp

0.000...0015e+xx の形式 atcoder.jp

どれも AC になってくれて、えらいなと思いました。\(2^{27}\) 文字出力しています*1。これを超えると OLE です。

話題 1

環境によっては、指数部が大きすぎると、無限精度で行った場合と異なる結果が返ってきそうです。 実数を入力として受け取り、double としてパースし、それを出力するプログラムを考えます。

#include <cstdio>
int main() {
  double x;
  scanf("%lf", &x);
  printf("%f\n", x);
}

これに対して、2000...00e-xx の形式で、2.0 と等しいはずの値を与えてみます。

% echo 20.0e-1 | ./a.out
2.000000

% (e=19999; printf 2; repeat $e printf 0; echo e-$e) | ./a.out
2.000000

% (e=20000; printf 2; repeat $e printf 0; echo e-$e) | ./a.out
20.000000

% (e=20306; printf 2; repeat $e printf 0; echo e-$e) | ./a.out
19999999999999999720621195205129155434005283676252727750499321471767131705345487698129692828457921333572758560785309230786706345700504206672551904741230794021461383329378750357138079702146292679283246532142253440022040339106608037192915625377123894402342976922345843644278133859702564244005353335500042141696.000000

% (e=20307; printf 2; repeat $e printf 0; echo e-$e) | ./a.out
inf

わーびっくり。AtCoder のコードテストの環境では再現しませんでした。 処理系によるのかな?

Because computers are finite, C++ implementations are inevitably limited in the size of the programs they can successfully process.

と規格にも書かれているし、そういうのも許されているのかもしれません? 上記は入力側じゃなくてプログラム側の話な気もするので、やっぱり許されてないかもしれません?

ソースコード中にそうしたリテラルを書いた場合はちゃんと出力してくれました。

% cat <<EOF | g++-10 -xc++ - -O9 -o fl && ./fl
cmdand heredoc> #include <cstdio>
cmdand heredoc> int main() {
cmdand heredoc>   double x = $(e=100000; printf 2; repeat $e printf 0; echo e-$e);
cmdand heredoc>   printf("%f\\n", x);
cmdand heredoc> }
cmdand heredoc> EOF
2.000000

cmdand heredoc> もプロンプトですが、syntax highlighting がこわれてしまった)

話題 2

ちなみに、浮動小数点数には 16 進のリテラルもあって、0xH.HHHp+D の形式です。H の部分を 16 進数で書き、D にはビットシフトの幅を書きます。 たとえば、\(2.75 = (2^3+2^1+2^0)/2^{-2}\) なので、2.75 の代わりに 0xbp-30.x1.6p+1 などと書けます。

printf("%a\n", x) などとすればこの形式で出力できますが、ジャッジはこの形式には対応していないようでした(かなしいね*2)。

atcoder.jp

おわり

今回は ABC 130 C のジャッジちゃんと遊びましたが、別のジャッジちゃんは別の挙動をするかもしれません。 誤差まわりについても遊んでみたいかもしれませんね。

*1:出力内容によっては 1 文字少ないかもです。

*2:別にかなしくはなくない?

bitset 高速化が定数倍高速化という風潮に対するお気持ち表明

「bitset で 64 倍高速化できます」「64 は定数なのでオーダーは変わりません」のような説明は、競プロ界隈でたびたび見かけます。

AGC の解説 とかも。

「定数倍高速化を軽視する風潮に対するお気持ち表明」や「競プロのような入力に上限がある場合でもオーダーを信仰する風潮に対するお気持ち表明」などは、一旦おいておきます。

今回のメインは、「長さ \(n\) の bit 列の bitwise AND/OR や shift にかかる時間は \(\Theta(n)\) か?」(あるいは、そう仮定するのは妥当か?)という話です。

補足

オーダー記法に慣れていない人が以下で怯えないように、さくっとまとめておきます。

  • \(\Theta(f(n))\):ちょうど \(f(n)\) のオーダー。
    • 競プロ界隈で言うところの \(O(f(n))\) に近そう。
  • \(O(f(n))\):高々 \(f(n)\) のオーダー。
    • \(O(n^3)\) と書いたとき、実際には \(\Theta(n^2)\) かもしれない。
  • \(\Omega(f(n))\):少なくとも \(f(n)\) のオーダー。
  • \(o(f(n))\):\(f(n)\) より真に小さいオーダー。
  • \(\omega(f(n))\):\(f(n)\) より真に大きいオーダー。

小文字の方は今回は使いません。

導入

次のような問題を考えてみます。

\(0\) または \(1\) からなる長さ \(n\) の列 \(a = (a_0, a_1, \dots, a_{n-1})\) が与えられます。これが回文か判定してください。

  1. for \(i \in \{0, 1, \dots, n-1\}\) do
    1. if \(a_i \ne a_{n-i-1}\) then
      1. return no
  2. return yes

愚直には、このようなアルゴリズムが考えられるわけですが、これの計算量はどうなるでしょう? 何の仮定も無しに \(O(n)\) 時間(+ worst \(\Theta(n)\) 時間)だと言い切っていいでしょうか? \(i\) のインクリメントや \(n-i-1\) の計算、\(a\) へのアクセスが \(O(1)\) 時間でできるものと仮定しているはずです(もっと言えば、1 bit の等価判定を \(O(1)\) 時間と仮定しているはずです。さすがにこれくらいは仮定したいですが)。

読者の中には、「\(a\) を順方向と逆方向にシーケンシャルアクセスするだけでいいからそれらの仮定はいらないだろう」と言う人もおられるかもしれない*1ので、別の例も出します。

\(0\) 以上 \(2^n\) 未満からなる長さ \(n\) の整数列 \(a = (a_0, a_1, \dots, a_{n-1})\) が与えられます。その後、「与えられた \(i\) と \(j\) に対して \(a_i+a_j\) の値はなんですか?」という形式の \(m\) 個の質問に答えてください。

これは \(O(n+m)\) 時間と言えてほしいですが、何の仮定も無しには言えません。 そもそも各整数をどういう形式で表すのかについても触れていません。 「ひとつの整数を表すのに \(2^n -1\) 個の連続する bit を用い、そのうち(連続するとは限らない)\(x\) 個の bit が \(1\) で残りが \(0\) のとき、\(x\) を表すとする」などと言われては大変扱いにくいので、もう少しマシなものを要求したいです。 コンピュータで扱いにくいものの話をしてもしょうがないので、2 進法で表すのがいいかなとなります*2

また、整数 \(i\) を用いて \(a_i\) に \(O(1)\) 時間でアクセスできることが要求されます。 \(a\) の長さは \(n\) であることから、\(\log(n)\) bits 程度の整数は \(O(1)\) 時間で扱えないとお話にならないことがわかります。 \(a_i+a_j\) の計算もしなきゃなので、\(\log(n)\) bits 程度の加算も \(O(1)\) 時間でしたいです。

word RAM の話

こうした要望に合った計算モデルとして、word RAM というのがあります。次のような仮定です。

  • \(U=2^u\) 個の word からなる配列をメモリと呼ぶ。
    • 入力はこのメモリに収まる必要がある。
  • word をポインタとして扱える。
    • word を添字として用いてメモリ(の全範囲)にアクセスできるということ。
  • 長さが \(O(1)\) の連続する word に対する読み書き、四則演算 (+, -, *, /, %)、ビット演算 (&, |, ^, !, <<, >>), 比較演算 (==, <, etc.) は \(O(1)\) 時間でできる。
    • ! はビット反転。C とかで言うところの ~

word によってメモリにアクセスできることから、word は最低でも \(u\) bit 必要です。 word の bit 長を \(w\) とすると \(w\ge u\) です。

これらから、連続する \(u\) bit に対する各種演算が \(O(1)\) 時間でできることが従います。 ただし、「メモリから飛び飛びの \(u\) 個の bit を集めてきて、結合してひとつの整数と見なし、それらを用いた演算を行う」ということは \(O(1)\) 時間でできるとは言っていないことに注意してください。

本題

さて、元々の話に戻ります。

長さ \(n\) の bit 列に対するビット演算の計算量は?

連続する \(n\) 個の bit 列、すなわち \(n/w\) word に対する演算となります。 \(w\) は最低でも \(\log(n/w)\) bit 必要なことから \(w \ge \log(n/w)\)、すなわち \(w\cdot 2^w \ge n\) となります。

式変形がんばるよ \(w\cdot 2^w = n\) の形の式を解くのは初等的にはつらいです。Lambert の \(W\) 関数 というのがあります。 一般には複素数で定義されますが、ここでは実数についてのみ考えます。 \[ ye^y = x \implies y = \begin{cases} W_0(x), & \text{if }x \ge 0; \\ W_{-1}(x), & \text{if } -1/e\le x\lt 0. \end{cases} \] ここでは \(x\ge 0\) のみ気にするので、\(W_0\) のみ考えます。 \[ \begin{aligned} w\cdot 2^w &\ge n \\ w\cdot 2^{w\cdot\log_2(e)\cdot\log_e(2)} &\ge n \\ w\cdot e^{w\cdot\log_e(2)} &\ge n \\ w\cdot\log_e(2)\cdot e^{w\cdot\log_e(2)} &\ge n\log_e(2) \\ w\cdot\log_e(2) &\ge W_0(n\log_e(2)) \\ w &\ge W_0(n\log_e(2))\cdot\log_2(e). \end{aligned} \]

また、\(W\) 関数について、以下のことが知られています*3。 \[ x\ge e \implies \ln(x)-\ln(\ln(x)) \le W_0(x). \] これらから、\(w\in \Omega(\log(n/\log(n)))\) が言えそうです。 がんばりおわり。

不等式を示すだけなら、\(W\) を使わなくてもできるかもです。

結局、\(w\in \Omega(\log(n/\log(n)))\) が言えたので、\(n/w \in O(n/\log(n/\log(n)))\) となります。 すなわち、bitset 高速化は(word RAM の文脈では)\(\log(n/\log(n))\) 倍の改善と言えそうです。

あるいは、入力が bit 列以外にもあるなどでメモリが \(\Omega(n)\) word 必要なら、\(w\in\Omega(\log(n))\) とできるので、\(\log(n)\) 倍の改善と言えそうです。

(追記:2021/06/11)\(O(n^2/w) = O(n^2/\log(n))\) になるやつとかを考えていました。メモリを \(\Theta(n)\) 使って \(O(n/w)\) だと、\(O(n)+O(n/w)=O(n)\) で得できなさそう?

おまけ

別のモデルの話

たとえば、もっとやばい計算モデルを仮定して、「任意のサイズの整数の四則演算は \(O(1)\) でできます」というのがある程度の妥当性をもって言えれば、bitset 高速化は \(n\) をひとつぶん落としたと主張できそうです。

word RAM の元々の論文*4では、

これこれの論文では長さ \(n\) の整数配列 \(a\) のソートを \(O(n)\) 時間で達成したと述べているが、\(n^2 \log(\max a)\) の長さ(ビット?)の除算を \(O(1)\) 時間でできると仮定した上でのアルゴリズムであり、やばいだろ

といった旨のことが書かれていました。

(追記)上の式では \(w\) の下限にしか触れていないので、別に \(w=n\) のようにワードサイズがやばいモデルを仮定することもできそうです。が、ちょっと乱暴だと思います。また、メモリを無駄に使うことで「\(w\) がもっと必要でしょう」と主張することもできそうですが、無駄に使ったメモリのぶんの時間のせいで損をしそうです。

www.cs.cmu.edu

今回の内容に関連する講義資料です。word RAM 以外にも、Turing machine論理回路を用いた場合に計算量はどうなるか?といったことが書かれています。

許容する演算の話

word RAM においては word size の四則演算を \(O(1)\) と言っていますが、乗算や除算はこれに含めないとする立場もあるようです。 これは、特定の条件を満たす回路によって加算は実現できるが乗算は実現できないという事情によるようです。 \(n\) bit の乗算が \(M(n)\) 時間でできれば、Newton 法によって \(n\) bit の除算も \(O(M(n))\) 時間でできるので、乗算を仮定するなら除算も仮定されそうです。

また、立っている最上位 bit の位置は、乗算を許せば次の手法で \(O(1)\) 時間で得られます。

rsk0315.hatenablog.com

これにより、立っている最下位 bit が n & -n で得られることと合わせて、立っている最下位 bit の位置も \(O(1)\) 時間で得られます。

word RAM 上のアルゴリズムの話

こういう記事があります。

qiita.com

まとめ

\(O(n/w)\) などの \(w\) は定数ではなくて、入力サイズ \(n\) に依っていると見なす立場もあるよ、という話でした。

マシンが入力に依ることに納得がいかない人は、「大きいサイズの問題を、メモリの少ない昔の PC で解くのは不可能」みたいなことを考えてみるといいかもしれません? 大きなメモリをちゃんと扱うことができるマシンであれば、それに伴っていくらか大きい値も扱えるようになる、といったような気持ちに近い気がします。

おわり

log は定数学派の人にとっては bitset 高速化は定数倍高速化ということがわかりました。いかがでしたか?

*1:「インクリメントはならし定数時間なので〜」とか思った人もいるかも?

*2:当然、文脈によってはコンピュータに縛られる必要はありませんが、ここでは一応競プロに近い文脈なので。

*3:Hoorfar, A., & Hassani, M. (2008). Inequalities on the Lambert W function and hyperpower function. J. Inequal. Pure and Appl. Math, 9(2), 5-9.

*4:Fredman, Michael L., and Dan E. Willard. "Surpassing the information theoretic bound with fusion trees." Journal of computer and system sciences 47.3 (1993): 424-436.

素数判定するよ

\(n\) 以下の素数の個数を求める関数 \(\pi(n)\) を持っていたとします。 このとき、\(n\) の素数判定は \(\pi(n)-\pi(n-1) = 1\) かどうかでできます。

fn is_prime(n: usize) -> bool {
    n > 0 && prime_pi(n) - prime_pi(n - 1) == 1
}

たとえば \(O(n^{3/4}/\log(n))\) でできます*1

お粗末さまでした。

おまけ

素数の数え上げが気になる人向けの記事:

rsk0315.hatenablog.com

*1:素直に \(O(n^{1/2})\) 時間の試し割り法をしましょう。

n 番目の素数を求めるよ

軽めの記事です。

\(n\) 番目の素数を \(p_n\)、\(n\) 以下の素数を \(\pi(n)\) と書くとします。 このとき、\(p_n\) は \(\pi(p_n) \ge n\) かつ \(\pi(p_n-1) \lt n\) を満たします。 よって、これは二分探索で求められます。

上限を決めるのが大変そうですが、指数探索っぽくすればいいです。 指数探索というのは、上限を \(2, 4, 8, \dots\) と倍々に増やしていき、true/false の境界を跨いだらその範囲で二分探索を行う方法です。 今回は、自明に \(p_n \gt n\) なので、上限の初期値はそれに基づいて決めてよいです(別に \(2\) から始めてもよいです)。

fn main() {
    assert_eq!(nth_prime(3), 5);
    assert_eq!(nth_prime(4118054813), 99999999977);
}

fn nth_prime(n: usize) -> usize {
    let mut lo = n;
    let mut hi = 2 * n;
    let small = |i| prime_pi(i) < n;
    while small(hi) {
        hi *= 2;
    }
    while hi - lo > 1 {
        let mid = lo + (hi - lo) / 2;
        if small(mid) {
            lo = mid;
        } else {
            hi = mid;
        }
    }
    hi
}

ここで、prime_pi別の記事 で書いたやつを使うと、\(O(n^{3/4}/\log(n))\) time です。 また、\(p_n = O(n\log(n))\) が知られているので、全体で \(O( (n\log(n))^{3/4}/\log(n\log(n))\cdot \log(n\log(n)) = O( (n\log(n))^{3/4})\) time です。

お粗末さまでした。

眠れない夜は素数の個数でも数えましょう

別に、おふとんから出たくない朝とかに数えてもいいと思います。

\(n\) 以下の素数の個数を数えたいときどうしますか? 各整数 \(2\le i\le n\) に対して素数判定して数えます? 素数判定を試し割り法でやると \(\Theta(n\sqrt{n})\)、線形篩*1を使えば \(\Theta(n)\) とかですね。 これを \(o(n)\) でできたらうれしくないですか?という話です。

今回は \(O(\sqrt{n})\) space, \(O(n^{3/4}/\log(n))\) time の話がメインです。 定数倍があまり軽くない \(O(n^{2/3})\) space, \(O(n^{2/3})\) time の話も少し触れます。

(追記:2021/06/05)少し難しめかもですが、アイディア自体は水色もあれば理解できると思います。

前提知識

以下、\(n\) 以下の個数の素数を \(\pi(n)\) と書きます。

導入

素数の個数を数えるにあたって、素数を実際に列挙する必要はないというのが重要です。 実際に列挙していては、「与えられた素数の次の素数はなに?」というのに \(O(1)\) 時間で答えられても、\(O(n/\log(n))\) 時間かかってしまいます*2

ここで、Eratosthenes の篩について考え直します。次のようなアルゴリズムです。

  1. \(1\) から \(n\) の bool 配列 \(a_i\) を用意する。\(a_1\) のみ false とする。
  2. \(i = 2, 3, \dots, \lfloor\sqrt{n}\rfloor\) に対して、次の処理をする。
    1. \(a_i\) が false なら次の \(i\) に移る。
    2. \(j = i, i+1, \dots, \lfloor n/i\rfloor\) に対して、\(a_{ij}\) を false にする。
let mut is_prime = vec![true; n + 1];
is_prime[1] = false;
for i in (2..=n).take_while(|&i| i * i <= n) {
    if !is_prime[i] { continue; }
    for j in i..=n / i {
        is_prime[i * j] = false;
    }
}

ここで、各 \(i\) に対して 2.2 のステップでいくつの \(a_{ij}\) が true から false になるかについて考えます。 この個数をうまく管理できれば、篩い終えた後にいくつ true になっているか、すなわち素数の個数がわかることになります。

この手法は、界隈では Lucy DP などと呼ばれることもあるようです*3

解説

(追記 2021/05/19)画像とかをつけるとわかりやすくなる気がしたので、そのうちつけると思います。

\(S(v, p)\) を、\(i=p\) まで処理を終えた後に \(a_j\) が true である \(2\le j\le v\) の個数と定義します。 \(p\cdot p \gt v\) であれば \(i=p-1\) のときから配列は変化しないので、\(S(v, p) = S(v, p-1)\) とわかります。\(p\) が素数でない場合も同じです。

\(p\) が \(p\cdot p\le v\) なる素数の場合はむずかしいのでちょっと考えます。 \(j = p, p+1, \dots, \lfloor v/p\rfloor\) に対して \(a_{pj}\) が false になりますが、\(j\) が \(p\) 未満の素数で割り切れるものはすでに篩われています。よって、篩われる個数は「\(p, p+1, \dots, \lfloor v/p\rfloor\) のうちで篩われていないもの」であり、これは \(S(\lfloor v/p\rfloor, p-1)-S(p-1, p-1)\) です。すなわち、次が成り立ちます。 \[S(v, p) = S(v, p-1) - (S(\lfloor v/p\rfloor, p-1) - \pi(p-1) ).\]

さて、\(\pi(n)\) は \(S(n, \lfloor\sqrt{n}\rfloor)\) と等しいので、これを求めたいわけですが、\(S(v, \bullet)\) を求めるためには(その時点で見つかっている素数の個数と)\(S(\lfloor v/\bullet\rfloor, \bullet)\) の形で表せるものさえわかっていれば十分とわかります。 \(\lfloor \lfloor n / p\rfloor / q\rfloor = \lfloor n / pq\rfloor\) が重要です*4

よって、\(v = 1, 2, \dots, \lfloor\sqrt{n}\rfloor\) と \(v = \lfloor n/1\rfloor, \lfloor n/2\rfloor, \dots, \lfloor n/\lfloor\sqrt{n}\rfloor\rfloor\) についてのみ管理します(サイズ \(O(\sqrt{n})\) の配列ふたつで管理できます)。 以下のコード中では \(\texttt{smalls}[i] = S(i, p)\)、\(\texttt{larges}[i] = S(\lfloor n/i\rfloor, p)\) とします。

fn main() {
    assert_eq!(prime_pi(10), 4);
    assert_eq!(prime_pi(100), 25);
    assert_eq!(prime_pi(100000000000), 4118054813);
}

fn prime_pi(n: usize) -> usize {
    let n_sqrt = floor_sqrt(n);

    let mut larges: Vec<_> = (0..=n_sqrt).collect();
    for i in 1..=n_sqrt { larges[i] = n / i - 1; }
    let mut smalls: Vec<_> =
        (0..n / n_sqrt).map(|x| x.saturating_sub(1)).collect();

    for p in 2..=n_sqrt {
        if p < smalls.len() {
            if smalls[p] <= smalls[p - 1] { continue; }
        } else {
            if larges[n / p] <= smalls[p - 1] { continue; }
        }
        let pc = smalls[p - 1];
        let q = p * p;
        for i in (1..=n_sqrt).take_while(|&i| n / i >= q) {
            // vi = n / i
            // dp[n / i] -= dp[n / i / p] - pc;
            let ip = i * p;
            let cur = *larges.get(ip).unwrap_or_else(|| &smalls[n / ip]) - pc;
            *larges.get_mut(i).unwrap_or_else(|| &mut smalls[n / i]) -= cur;
        }
        for i in (1..smalls.len()).rev().take_while(|&i| i >= q) {
            // vi = i
            // dp[i] -= dp[i / p] - pc;
            let cur = smalls[i / p] - pc;
            smalls[i] -= cur;
        }
    }
    larges[1]
}

fn floor_sqrt(n: usize) -> usize {
    if n <= 1 {
        return n;
    }
    let mut lo = 1;
    let mut hi = n;
    while hi - lo > 1 {
        let mid = lo + (hi - lo) / 2;
        match mid.overflowing_mul(mid) {
            (x, false) if x <= n => lo = mid,
            _ => hi = mid,
        }
    }
    lo
}

計算量解析

i に関するループがふたつありますが、ひとつめの計算量は次の値で抑えられます。以下の \(i\) は p に対応します。 \[ \begin{aligned} \int_2^{\sqrt{n}} \min\{\sqrt{n}, n/i^2\}\, di &= \int_2^{\sqrt[4]{n}} \sqrt{n}\,di + \int_{\sqrt[4]{n}}^{\sqrt{n}} n/i^2\,di\\ &= \sqrt{n}\cdot\sqrt[4]{n} + O(1) + (\sqrt[4]{n}-1)\cdot\sqrt{n} \\ &= O(n^{3/4}). \end{aligned} \]

ふたつめも。 \[ \begin{aligned} \int_2^{\sqrt{n}} \max\{n-i^2, 0\}\, di &= \int_2^{\sqrt[4]{n}} (\sqrt{n}-i^2)\, di + \int_{\sqrt[4]{n}}^{\sqrt{n}} 0\,di \\ &= 1/3\cdot (2n^{3/4} - 3\sqrt{n} + 1) \\ &= O(n^{3/4}). \end{aligned} \]

合わせて \(O(n^{3/4})\) です。実際には、素数の分布から \(O(n^{3/4}/\log(n))\) になっている気がします。

工夫

オーダーの改善ではないはずで、定数倍高速化だと思います。

上記のアルゴリズムで、篩う \(i\) を \(\lfloor\sqrt{n}\rfloor\) までではなく \(\lfloor \sqrt[4]{n}\rfloor\) までにしてみます。 このとき、篩われずに残っている合成数は、たかだか 3 つの素数の積で表せることがわかります*5。4 つかけると \(n\) を越えちゃうので。 このことから、全体をふたつの処理に分けて行うことが考えられます。

  • \(\lfloor\sqrt[4]{n}\rfloor\) 以下の素数で篩い、篩われた個数を引く。
  • たかだか 3 つの素数の積で表せる合成数の個数を引く。

合成数の個数を数える方について考えます。 篩われずに残っている整数(以前篩うのに使ったのは除く)を昇順に見ていきます。\(i_1, i_2, \dots\) とします。 各 \(k\) に対して、\(i_k\) の倍数である \(n\) 以下の整数の個数を足します。 また、\(k'\lt k\) に対して \(i_k\cdot i_{k'}\) の倍数である整数はすでに数えられているので、引きます。

ここで、\(i_k\) の倍数である \(n\) 以下の整数は、\(S(\lfloor n/i_k\rfloor, \lfloor\sqrt[4]{n}\rfloor) - \pi(i_k -1)\) と表せます。\(i_k i_{k'}\) の方は \(S(\lfloor n/i_k i_{k'} \rfloor, \lfloor\sqrt[4]{n}\rfloor) - \pi(i_{k'} -1)\) です。

実装

アイディアは上記の通りですが、実装上の工夫が少しあるので下で説明します。

変数の概要は次の通りです。

  • \(\texttt{smalls}[i] = S(i, p)\)
  • \(\texttt{roughs}[i]\):消されていない \(i\) 番目の整数
    • \(1\) は残る、篩うのに使った素数は消える
    • サイズは \(\texttt{s}\) で管理される
  • \(\texttt{larges}[i] = S(\lfloor n/\texttt{roughs}[i]\rfloor, p)\)
  • \(\texttt{pc} = \pi(p-1)\)
  • \(\texttt{skip}[i]\):\(i\) が篩われていれば true
    • ただし偶数については更新しない

\(p\) の変化に伴って更新されるので場所によって ±1 がありますが、各ループの最後で上記の定義に合うように管理します。 最初のループの開始前においては、\(p = 2\) での値です。

fn main() {
    assert_eq!(prime_pi_fast(100000000000), 4118054813);
}

fn prime_pi_fast(n: usize) -> usize {
    if n <= 3 {
        return n.saturating_sub(1);
    }
    let v = floor_sqrt(n);
    let mut smalls: Vec<_> = (0..=v).map(|i| (i + 1) / 2).collect();
    let mut s = (v + 1) / 2;
    let mut roughs: Vec<_> = (0..s).map(|i| 2 * i + 1).collect();
    let mut larges: Vec<_> =
        (0..s).map(|i| (n / (2 * i + 1) + 1) / 2).collect();
    let mut skip = vec![false; v + 1];

    let mut pc = 0;
    for p in (3..=v).step_by(2) {
        if skip[p] { continue; }
        let q = p * p;
        pc += 1;
        if q * q > n { break; }
        skip[p] = true;
        for i in (q..=v).step_by(2 * p) {
            skip[i] = true;
        }
        let mut ns = 0;
        for k in 0..s {
            let i = roughs[k];
            if skip[i] { continue; }
            let d = i * p;
            let x = if d <= v {
                larges[smalls[d] - pc]
            } else {
                smalls[n / d]
            };
            larges[ns] = larges[k] + pc - x;
            roughs[ns] = i;
            ns += 1;
        }
        s = ns;
        let mut i = v;
        for j in (p..=v / p).rev() {
            let c = smalls[j] - pc;
            let e = j * p;
            while i >= e {
                smalls[i] -= c;
                i -= 1;
            }
        }
    }

    let roughs = roughs;
    let pc = pc;

    let mut res: usize = larges[0] + (s + 2 * (pc - 1)) * (s - 1) / 2
        - larges[1..s].iter().sum::<usize>();

    for l in 1..s {
        let q = roughs[l];
        let m = n / q;
        let e = smalls[m / q] - pc;
        if e <= l { break; }
        let t: usize = roughs[l + 1..=e].iter().map(|&r| smalls[m / r]).sum();
        res += t - (e - l) * (pc + l - 1);
    }
    res
}

fn floor_sqrt(n: usize) -> usize { /* 上と同じなので省略 */ }

実装の説明

冷静になると自明かも。わからなければ、下記の説明を読みつつ、各箇所で値を適宜出力したりしながら確かめるとよいと思います。

roughs[ns] の更新

roughs[..s] が篩われていない整数です。 ns が次のループにおける s です。 ns = 0 で初期化しておき、roughs[ns] に値が入ったら ns += 1 していきます。前から詰め直しているということですね。

larges[smalls[d] - pc] についての説明

smalls[d] - pc は「篩われずに残っている整数のうち d は何番目か?」を意味します*6

smalls[d] - pc は「d 以下にある篩われていない整数」から「篩うのに使った整数」を除いたものの個数です。 ただし、この段階ではまだ pc は \(\pi(p -1)-1\) であることに注意です。これに、\(1\) のぶんの個数を足して、所望の値であることがわかります。 roughs[それ] が篩われていない d 番目の整数であることと larges の定義から、結局 \(n/d\) 以下で篩われていない整数の個数になります。

smalls[i] の更新

j = i / p を保って更新されるようにします。 i >= j * p である間 smalls[i] -= smalls[j] - pc で更新され、smalls[i] -= smalls[i / p] - pc になっています。

res の初期化

初期化の際に、(上で説明した)\(i_k\) の倍数を先に引いています。等差数列の公式などを思い出すと吉です。

res の更新

  1. q (\(i_{k'}\)) を固定
  2. q より大きい r (\(i_k\)) を固定
  3. n / (q * r) として表せる篩われていない整数の個数を求める

\(\sqrt[4]{n} \lt i_{k'} \lt i_k\) から、\(\lfloor n / i_k i_{k'}\rfloor \le \lfloor\sqrt{n}\rfloor\) であり、m / r = n / (q * r)smalls にアクセスしても問題ありません。

3. で求めた個数から、(前半パートですでに数えた)\(\sqrt[4]{n}\) 以下の素数を含む個数を除く処理が -= (e - l) * (pc + l - 1) です。 \(\lfloor n / i_k i_{k'}\rfloor\) として表せるうちで数えたい最小値は \(i_k'\) なので、2. で固定するのは \(\lfloor n / i_{k'}^2\rfloor\) 以下の範囲です。

計算量解析

前半パートでは、各 \(i \le n^{1/4}\) の素数に対して \(O(\sqrt{n})\) 程度の処理をするので \(O(n^{3/4}/\log(n))\) time です。篩の部分も \(O(\sqrt{n}\log(\log(n))) \subset o(n^{3/4}/\log(n))\) です。 後半パートでは、雑に見積もって \(n^{1/4} \lt i \le n^{1/2}\) の各素数に対して \(O(n/i^2)\) 程度の処理をしています。 \[ \int_{n^{1/4}}^{n^{1/2}} \frac{n}{i^2}\, di = n^{3/4} - n^{1/2}. \] であり、素数の分布から \(O(n^{3/4}/\log(n))\) になりそうです。

\(n = 10^{13}\) でも 1 秒前後くらいで終わりました。 \(\log(n)\) で割れる部分のちゃんとした(お気持ちでない)証明ができたら追記します。

その他

\(\sqrt[3]{n}\) とかで区切って \(O(n^{2/3})\) 時間にできるのかな? 区切り方を変えれば \(O(n^{3/4-\varepsilon})\) にはできそうな気もするけど、定数倍とか諸々のことはわかりません。

Dirichlet 積などを使った方法 とかもあります。軽めの定数倍では実装できませんでした。

参考文献

おわり

おやすみなさいにゃ。

(記事を書いていたら眠れなくなりました)

*1:実装はこれとか素因数分解の一意性をうまく使いつつ配る DP をします。

*2:素数定理から。

*3:Project Euler での Lucy_Hedgehog 氏の投稿 に由来するようです。

*4:ABC 141 D の解説 とかに証明があります。

*5:\(\lfloor\sqrt[3]{n}\rfloor\) までにすればたかだか 2 つにできます。

*6:rank/select 演算を知っている人は、それに似たことを考えるとよいと思います。