えびちゃんの日記

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

素数の数え上げのスライドについて

最近スライドを公開しました:

動機について

やや前に Beamer でスライドを作る機会があったのと、最近 LuaLaTeX を触ってみたというのがあって、LuaLaTeX + Beamer でスライドを書きたいなという欲がありました。

素数の数え上げや乗法的関数の和は、一部界隈では知られていそうだったものの、あまり浸透していない印象があったので勉強したいと思っていました。 そこで、それをスライドに書くことで一石たくさん鳥にしようと思いました。文字ベースの記事だと説明がごちゃつきそうに感じたというのもあります*1

スライドに関して

スライドのこだわりポイントや、作っていてうれしかった点などをぽんぽん挙げていきます(淡々と空行区切りで書いたところ思ったより読みにくくなった感があります)。

暗めの赤と青の配色は、大学時代に学生たちから畏れられていた教授がよく使っていたテーマを真似して作ったものです*2

Beamer に詳しい人は知っているかもしれませんが、スライドの以下のような部分がリンクになっていて、ページ遷移が容易なのがうれしいです。

  • ヘッダ部にあるセクション名
  • ヘッダ部にある〇
  • フッタ部左側のタイトル

同じトピックのスライドに関して、タイトルに「I, II, III, IV, ...」と振っている箇所が複数ありますが、これはそういう制御綴を置いておくことで自動で採番してくれるのでうれしいです。

\newcounter{slidetopic}
\renewcommand{\theslidetopic}{\stepcounter{slidetopic}\Roman{slidetopic}}

\setcounter{slidetopic}{0}
\begin{frame}
  \frametitle{タイトル \slidetopic}  % タイトル I
\end{frame}
\begin{frame}
  \frametitle{タイトル \slidetopic}  % タイトル II
\end{frame}
% ...

右下にあるページ番号も(当然)自動採番ですが、\appendix と書いた後は appendix 用のものと切り替えられるらしいのを知ってうれしくなりました。

\addtobeamertemplate{navigation symbols}{}{%
  \usebeamerfont{footline}%
  \usebeamercolor[fg]{footline}%
  \hspace{1em}%
  \makeatletter%
  \ifbeamer@inappendix{%
    \textsc{appendix} {\insertframenumberinappendix} / {\insertappendixframenumber}%
  }\else{%
    {\insertframenumber} / {\insertmainframenumber}%
  }\fi%
  \makeatother%
  \hspace{0.8em}%
}

小文字の appendix だと配置が崩れたので small caps を使いました。

\(+=\) ではなく \(\xleftarrow{+}\) と書くのは、\(\gets\) で代入を意味する擬似コードではいい感じの見栄えなんじゃないかなぁと自画自讃しています。 そこはいいと思うんですが、「記法について」のスライドがちょっと突飛かもと思いました。「記法と言えば prelim で、prelim と言えば最初でしょ」みたいなノリだったと思います。

木の部分は気に入っていて、TikZ はえらいなあと思いました。

\usetikzlibrary{
  graphs,
  graphdrawing,
}
\usegdlibrary{trees}

% ... in {tikzpicture}
\graph[tree layout, sibling sep=0pt]{
  1 -> {
    2 -> {
      4 -> {
        8 -> 16,
        12,
        20
      },
      6 -> 18,
      10,
      14
    },
    3 -> {
      9,
      15
    },
    5,
    7,
    e/...,
    19
  }
};

辺ラベルは、白い太文字を書いた上に黒文字をかぶせることで実現しています。

\usepackage{pdfrender}
\newcommand{\tpr}[1]{%
  \textpdfrender{
    TextRenderingMode=FillStrokeClip,
    LineWidth=3pt,
    FillColor=black,
    StrokeColor=white,
    MiterLimit=1
  }{#1}%
}

% ... in {tikzpicture}
\begin{scope}[every node/.style={font=\tiny}]
  \path
  (1) -- node {\tpr{2}} (2)
  (1) -- node {2} (2)
  (1) -- node {\tpr{3}} (3)
  (1) -- node {3} (3)
  % 各辺に対して同じようなことを書いた。Lua を使えばもっと楽にできたと思う。
  ;
\end{scope}

この構成方法(自身を最小素因数で割ったものが親になる)で作られた木に名前がついていないのか気になります。ないのかなあ。

冒頭の篩や付録ページの DP の差分スライドたちは Lua で実際に DP しながら値を計算して、その結果を出力しています*3。 こういうのが柔軟にできるのはうれしいなあという気持ちになります。Lua なのがうれしいかというと怪しいですが、十分ありがたいですね。

付録スライドに載せるようなコンテンツ(擬似コードなど)は、発表後にゆっくり読んでもらう側面が強いと思ったので、付録の擬似コードはめちゃくちゃに縮小してみました。 最初は、

function lucy()
|  foo
|  bar
+ ...
function lucy()
| ...
| baz
| qux
+ ...
function lucy()
| ...
| quux
| corge
+---

などのように分割することも考えたのですが、シンプルに読みにくいし書くのも大変だなとなってやめました。 スライドを PNG などの画像に変換して読んでいる人がいたらかわいそうなことになりますが、まぁいいかなと思いました。

改行位置には結構気をつけたつもりです。

おまけスライドにちょっとした問題を追記したのですが、「読者への課題です」と言って説明を放棄して迷宮入りするのは好ましくないと思ったので、文字を逆向きにしてちょっと読みにくくして載せる手法を使いました。小さい子向けのなぞなぞの本とかでこういうのがあったような気がしていますが、参考書とかではあまり見ないような気がします。 えびちゃんのセグ木スライドでも同様の手法を使った気がします。

スライドを更新した際に(スライドの TeX ファイルなどは private repository で管理しているのですが)commit hash をつけておくと何らかの捗りがある気がしたので、更新日時とともに載せてみました。LuaLaTeX で比較的簡単に実現できたのでうれしかったです。

ねこねこ勉強ぱーてぃというのは架空の勉強会の名前です。とりあえず仮置きで名前をつけていたところ、考え直すのを忘れて投げてしまいました。

承認欲求の満たされについて

うれしい。競プロ界隈でスライドを作っている人が少ししかいないかもという気持ちも少しあります。と見せかけてそんなことないかも?

うれしい。書いてみた甲斐があります。

その他

実験コードを書いて、計測をした後、「念のため擬似コードにしておくか」と思ってコードを読んでいたところ、実装に誤りに気づいてつらくなりました。\(p_{\pi} \gt n^{1/6}\) で判定するところを \(\pi \gt n^{1/6}\) と書いていました。たぶん log 倍悪くなっていそうですが、修正したらちゃんと速くなってくれたのでうれしかったです。

擬似コードで書く利点ってなんなんでしょう。Rust のコードをそのまま貼る方がうれしいかもとも思ったりはしました。 \(\text{\textbf{foreach }}v \in S\) とか \(x \xleftarrow{+} \sum_i f(i)\) のように数式混じりの表現をそのまま使っても違和感がなさげなのがちょっとうれしいかもと思ったりはしました*4。 特定の言語依存の表現だと、たとえば log(x, y) と書いたときどちらが底なのかわかりにくいとか、そういう嫌な点があるかなとも思いました。

どの言語にも依存していないが、どの言語に移植するにも手間がかかってしまうというのはありますね*5

擬似コードオブジェクト指向っぽく \(v.\text{push}(x)\) とか書くのが邪道なんじゃないかと思ったりしつつ、どう書くべきかわからなかったりしました。

あと \(\text{\textbf{for each}}\) のように 2 語にする方が気持ちいいんじゃないかという気持ちもあります。どうしよう。

\(\Theta(n^{2/3})\) の方は \(\Theta(n^{3/4}/\log(n))\) のものよりも定数倍が重くてしんどいかもという先入観があったのですが、意外と勝っていてうれしかったです。\(\Theta(n^{2/3}/\log(n)^{1/3})\) の方も速くてご満悦です。

最初は乗法的関数の \(\Theta(n^{2/3})\) の方は難しいと思ったので 2 回シリーズにしようかと思ったのですが、分かれちゃうのもアレかと思って、詰め込んでしまいました。

きもち

AtCoder でこの手のアルゴリズムが出題されるされないとかそういうのは抜きにして話します。 内容自体はそこまでレートが高くなくても読める内容だと思うので、読んでみてもらえるとうれしいかもしれないと思いつつ、書いた後だからそう思うだけで難しいのかもと思ったり、ぐるぐるしています*6

参考にしていたブログは最初は行間たっぷりに見えていたのですが、今では簡単に見えてしまったりしているので、えびちゃんの目はあてにならないです。

省略する。

とブログで書かれていた部分も、今では「そりゃ省略したくもなるよなぁ」と思ってしまうまであります。

なつかし

floor sum の解説記事を書いたときや、ACL math の解説記事を書いたときにあれこれ調べたのを思い出しました。数学のことを調べるのはやっぱりすきなんだろうなぁと思ったりしました。

おわり

おわりにゃ。次はなにの勉強会をしようかにゃぁ。

*1:実際、前に書いた記事はそういう感じになった気がします。

*2:Thank you! スライドも同様です。

*3:SATySFi でもそういうことができそうな感はありますね。

*4:下手に書くと計算量のごまかしや行間になりそうですが、それについてはスライド中で説明をしているつもりです。

*5:擬似コードってラテン語っぽい雰囲気があったりするでしょうか。

*6:あまり難しくないですよーと言って読んでもらおうとするのは、読めなかったときに落ち込んでしまいがちなので、あまりよくない気がしています。

不必要に浮動小数点数を使うの怖い

競プロにおいて、「何かしらの整数を大小比較して、左辺が大きければどうのこうの」みたいなことをしたい局面はそれなりにあります。 特に、境界付近でも正確に評価できる必要がある状況も多いです。 そうした状況で浮動小数点数を使うのはできるだけ避けたいですという気持ちがあります。

導入

まず、浮動小数点数は誤差が出ます*1

初心者のうちは「浮動小数点数では \(1\div 10\) を正確には計算できない」と聞いても、「たとえば

assert_eq!(1.0 / 10.0, 0.1);

とかはエラーにならないし、問題なく計算できているのでは?」みたいな誤解をすることもあります。 実際には、0.1 と書いた部分も計算前の時点で 0.1000000000000000055511151231257827021181583404541015625 のような値になっており、所望の状態にはなっていません。 1.0 / 10.0 もそういう値になっていたためエラーにならなかっただけですね。

  • そもそも表せない実数が(たくさん)ある
  • 表せる数同士の演算結果でも表せない数になることも(よく)ある

というのをちゃんと意識する必要があるかもしれません。 もちろん、表せなかった数同士で計算した結果がたまたま無限精度での計算結果と同じになることもあります。

assert_eq!(0.1 * 10.0, 1.0);

実際には、0.1000000000000000055511151231257827021181583404541015625 * 10.0 を計算しようとして、\(1.000000000000000055511151231257827021181583404541015625\) を表せずに 1.0 に丸められた結果、見かけ上は等しくなっているわけですね。

というわけで 1.0 / 10.0 * 10.0 == 1.0true ですが、たとえば 1.0 / 49.0 * 49.0 == 1.0false になります。 左辺は 0.99999999999999988897769753748434595763683319091796875 のような値になっています。

また、他によくある話として 0.1 + 0.2 != 0.3 というのもあります。 誤差をどう丸めるかにはいくつかのモードがありますが、「表せる数のうち最も近いものに丸める。ちょうど中間にある場合は、最も近い偶数の方向に丸める」という(よくある)ものを仮定すると、

0.1 => 0.1000000000000000055511151231257827021181583404541015625
0.2 => 0.200000000000000011102230246251565404236316680908203125
0.3 => 0.299999999999999988897769753748434595763683319091796875

となりますが、加減乗除については「無限精度での演算をしてから丸めた値と等しくなる」というように規定されているので、0.1 + 0.2、すなわち 0.100000000000000005551... + 0.200000000000000011102... を無限精度で計算してから丸めた 0.3000000000000000444089... のようになり、0.3 を丸めた 0.299999999999999988897... とは等しくなりません。「0.1 + 0.2 と書いたんだから、2 進法とかにはせずに書いた時点の値を無限精度で計算してから丸めてくれ」というのはきつそうです*2

0.1 + 0.2 => 0.3000000000000000444089209850062616169452667236328125

型について

doublef64 と言われる型でどのように数を表すかについて触れます。 IEEE 754 という規格で binary64 として定義されているものです*3。 実数の他に \(\infty\) に対応する値 Infinity や、非数(定義域外の値を関数に入力したときの返り値になったりする)に対応する値 NaN などもこれらの型には含まれますが、一旦それらは無視します*4

例として 3.141592653589793115997963468544185161590576171875 を挙げます。これは以下のようなビット列で表されます。

0 10000000000 1001001000011111101101010100010001000010110100011000

三つに分けて書いた部分をそれぞれ \( (s, e, b_{51} b_{50} \dots b_1 b_0)\) とおくと、次の式で表される値に対応します。 \[ (-1)^s \cdot 2^{e-1023} \cdot (1.b_{51} b_{50} \dots b_1 b_0)_{(2)}. \] つまり、以下のような感じです。 \[ \begin{aligned} {} &{\phantom{{}={}}} 2^{1024-1023} \cdot 1.1001001000011111101101010100010001000010110100011000_{(2)} \\ &= 2\cdot 1.1001001000011111101101010100010001000010110100011000_{(2)} \\ &= 2\cdot (1 + 2^{-1} + 2^{-4} + 2^{-7} + 2^{-12} + \dots + 2^{-48} + 2^{-49}) \\ &= 3.141592653589793115997963468544185161590576171875 \end{aligned} \]

というわけで、一番上の \(1\) がある桁から見て \(52\) 桁以上離れた \(1\) があるような数 (\(2^k\cdot 1.\underbrace{00\dots 0}_{52}1_{(2)}\)) や、十分大きい数 (\(2^{1024}\)) などは表せないことがわかります。このことから当然、有限小数でない実数も表せないことがわかります。\(2\) 進法において \(0.1_{(10)}\) は有限小数ではないため、\(0.1\) も表せません。

bartaz.github.io

こういうサイトがあって、オサレな UI で遊べます。

有限小数について

\(b\) 進法で既約分数 \(p/q\) が有限小数となる条件について考えておきます。簡単のため \(0\lt p/q\lt 1\) とします。

ある有限の \(k\) について \[ p/q = (0.x_1 x_2\dots x_k)_{(b)} \] とすると、両辺を \(b^k\) 倍して \[ b^k\cdot p/q = (x_1x_2\dots x_k)_{(b)} \] です。すなわち、\(b^k\cdot p/q\) が整数になる必要があります。\(p\) と \(q\) は互いに素なので \(b^k/q\) が整数になる必要があります。 このとき、(\(b\) と \(q\) の各素因数について約分していくみたいなことを考えると)\(q\) の素因数を \(b\) がすべて持っている必要があります。 逆に、\(q\) の素因数を \(b\) がすべて持っていれば十分です。

\(0.1 = 1/10\) に関しては、\(10\) の素因数 \(5\) を \(2\) が持っていないため、\(2\) 進法では表せないわけです。

嫌な例

単純な分数の比較

「\(x / y \le z / w\) か? (\(x, y, z, w\in \Z\cap[-10^2, 10^2]\), \(0\notin\{y, w\}\))」というのを知りたいケースはまぁよくあります。 これは、分母を払って \(xw\le yz\) とすれば整数で比較できます*5が、浮動小数点数で計算してしまう初心者は多いです。

ざーーっと調べてみた感じ、hack できるケースは見つからなかったのですが、コンテスト本番で不安になりながら書きたくはないので、整数でやりたいです。

大きい整数の比較

「\(xy \ge z\) か? (\(x, y, z\in\Z\cap[0, 10^{18}]\))」というのを知りたいケースはそこそこあるでしょう。x * y がオーバーフローするので整数でやりたくないと思ってしまう人が多そうです。

fn f(x: i64, y: i64, z: i64) -> bool {
    let x = x as f64;
    let y = y as f64;
    let z = z as f64;
    x * y >= z
}

たとえば、\(10^{18}-1 = 27027027\times 37000000037\not\ge 10^{18}\) ですが、f(27_027_027, 37_000_000_037, 10_i64.pow(18))true になってしまいます。 左辺が誤差で 1_000_000_000_000_000_000.0 になってしまい、等しくなってしまうんですね。

参考袋 を見つつ、\(x \ge \ceil{z/y}\) と変形します。\(y = 0\) に気をつけて、

fn f(x: i64, y: i64, z: i64) -> bool {
    if y == 0 { z == 0 } else { x >= (z + y - 1) / y }
}

のようになる気がします。

もっとも、i128 のような大きい整数型を使うという手もありますが、制約が \(10^{36}\) とかだったら?と考えると趣旨は伝わると思います。

小数点の移動

「小数点以下高々 \(9\) 桁の実数 \(0\lt x\lt 10^4\) が与えられます。\(10^9\cdot x\) を出力してね」というの を考えます。

適当に 1.0e9 倍して少し足しておけばいいでしょうみたいな気持ちになる人も多いと思うのですが、

fn f(x: f64) -> i64 { (1.0e9 * x + 1.0e-6).floor() as i64 }

などと書くとうまくいきません。"32.000000001" が与えられたとき、これをそのまま f64 に変換すると 32.00000000099999653002669219858944416046142578125 となり、32000000000.999996185302734375 となり、32000000000 となります。

文字列としてがんばって処理するか、.round() as i64 のようにするのが無難そうです*6

平方根の比較

これは有名。「\(\sqrt{a}+\sqrt{b}\lt\sqrt{c}\) ですか?」というやつです。

解説 を読んだり、maspy さんのツイート を読んだり。

累乗関連

Rust だと整数にも .pow() が用意されているのであんまりなさそうですが、C++ だと浮動小数点数にしか用意されていないせいか、整数の累乗をしたいときに浮動小数点数を介してしまう人が多いです。

「集合 \(S, T\) が与えられます。\(|S\cap T|\) は? (\(S, T\subseteq\{0, 1, \dots, 62\}\))」というのがあったとき、

// C++
long long x = 0, y = 0;
for (auto si: s) x += pow(2, si);
for (auto ti: t) y += pow(2, ti);
return __builtin_popcountll(x & y);

みたいなことを書いてしまう人もいます。powdouble を返すので、以下のような変換が暗黙になされています。

// x += pow(2, si);
x = (long long)((double)x + pow(2, si));

その結果、たとえば \(S = \{0, 53\}\) のようなとき、期待する値が \(9007199254740993\) なのに対し、誤差の関係で x == 9007199254740992 になってしまいます。 x |= 1LL << si のように書くのが無難な気がします。Rust においては x |= 1 << si と書いても問題ないです。

Fibonacci 数

\(n\) 番目 (\(1\le n\le 10^{18}\)) の Fibonacci 数を \(998244353\) で割ったあまりを求めてね、と言われて、\(F_n = (\sqrt{5})^{-1}\cdot (\phi^n-(-\phi^{-1})^n)\) を浮動小数点数演算で求めようとする初心者もいる気がします。

たとえば \(F_{72} = 498454011879264\) ですが、上式をそのまま計算すると、498454011879265.2 とかになってしまいそうです。

こっち はギャグの記事。

もっと単純な比較で hack できそうな例を出して初心者を脅かしたかったんですが、いいのが思いつきませんでしたね。

おまけ

これは「整数を使えばいいですよ」系とは違って、「こういう例もあります。嫌ですね」みたいなのを紹介するコーナーです。

誤差のせいで定義域を超えるやつ

敢えて作ったような例ですが、「\(x\) が与えられます。\(\sqrt{\sin(\pi x)}\) を出力してね」という問題を考えます。\(x = 2\) のときなどにうまくいきません。

use std::f64::consts::PI;

fn main() {
    println!("{:.50}", (2.0 * PI));
    println!("{}", (2.0 * PI).sin());
    println!("{}", (2.0 * PI).sin().sqrt());
}

以下のような出力になります。

6.28318530717958623199592693708837032318115234375000
-0.00000000000000024492935982947064
NaN

0.0 になってくれるでしょというのは思い込みなわけですね。

\(2\pi = \tau\) を表す定数 TAU も用意されていますが、TAU.sin().sqrt()NaN になります。

他にも、asin (\(\arcsin\)) の引数が 1.0 よりわずかに大きかったときも NaN になったりします。 定義域が \( (-\infty, \infty)\) でない関数*7を使うときは特に気をつけましょう。

意図しない感じに丸められるやつ

これも敢えて作った例ですが「\(1\le x\le 9007199254740995\) が与えられます。\(1/(9007199254740996-x)\) を出力してね」という問題を考えます。 x = 9007199254740995.0 とすると、x == 9007199254740996.0 となるように丸められ、1.0 / 0.0 == Infinity となります。

細かいやつ

\(1 / 3\) は 2 進法で有限小数にはならないので、\(\sqrt[3]{x} = x^{1/3}\) と思って計算すると少しずれることがあります。

println!("{:.050}", 5.0_f64.cbrt());
// 1.70997594667669705614798658643849194049835205078125

println!("{:.050}", 5.0_f64.powf(1.0 / 3.0));
// 1.70997594667669683410338166140718385577201843261719

競プロの範囲でこれが問題になることは少ないかもですが、頭の片隅に置いておいてもよいでしょう。

丸めに関して

加減乗除については、無限精度の演算結果を丸めたかのように振る舞ってくれます。

FMA (fused multiply-add; \(xy+z\)) というのがあって、丸めを一度しか行わないものとして定義されていますが、正しくない結果を返す実装も あるらしい です。

それ以外の関数 (\(\sin(x)\), \(\sqrt{x}\), \(x/y\), etc.) についてはどうでしょうか。 IEEE 754-2019 では required operations と recommended operations という枠が分かれていて、加減乗除などは前者、\(\sin\) などは後者に含まれています*8。 後者も無限精度の演算結果を丸めたかのようにするべきと書かれていますが、required ではないのでゆるっとした感じがする気がします。実際、GCC は一部を除いて特にちゃんとすることは目指していない らしいです

table maker’s dilemma とかを調べるとよいかも。Lindemann の定理 とかも関連していそう。 関連してそうな資料、これ とか これ とか。

詳しい人のツイート:

その他

もしかすると x < y + eps (eps = 1e-8) みたいにするといいとかを聞いたことがある人もいるかもしれません。 それはどうしても浮動小数点数を使う必要がある場面で適切に使うべきであって、整数に帰着できる問題ではやらない方が無難そうです。

根本的に間違った解法の場合でも「eps を変えれば通るかも」などと言って無駄に誤答を出すことにもつながりがちです。

また、幾何や幾何以外の問題などで分数を扱う場合においても、適当な素数 \(p\) を用いて、\(x/y\) を \(x\cdot y^{-1}\bmod p\) として表すとよいケースもあります。(問題 / 解説) (問題 / 解説)

こと浮動小数点数においては、初心者のうちに「なんかこれでうまくいきそう(いくつかのケースではうまくいった)」と思ったものには、反例がたくさんある気がします。 そうした方法を使って祈るよりは、(十分に勉強して祈らず使うか)整数に帰着する方法を考えるのが賢明そうです。

関連問題

浮動小数点数を避けてやりたいな〜となる問題を挙げておきます。浮動小数点数でやって問題があるかどうかは気にしていないので、問題がない場合でも整数でやる練習としてやるとよいでしょう。

他にもいい感じのがあったら募集しています。

計算に関して

docs.rs

浮動小数点数の数学関数の実装が載っています。 ハードウェアのサポートを得られない状況を想定して書かれていそうで、多くの状況ではハードウェア命令で処理されていそうな気がします。

近似とかをがんばったりしていそうです。

この辺は「不必要に〜」の文脈とは違いそうですが、ついでに載せちゃいました。

便利ちゃん?

誤差が減るように式変形してくれるらしいです。

herbie.uwplse.org

おわり

おわりです。

*1:誤差以外もいろいろ出そうと思えば出るかもしれません。

*2:コンパイル時の最適化によっては起きうる? オプションとかにもよるかも。C/C++ ではコンパイル時と実行時で精度や丸め方が変わるのを許容するみたいな話がどこかにあったような気がするのでまた調べます。最適化で誤差関連の挙動を変えてもいいよーみたいなオプションはあったはず。

*3:Rust はそうっぽいけど C++ は必ずしもそうではないかも。GCC ではサポートしていそう。そうじゃない処理系があるかは知らないけど規格上は許容されていそう。

*4:さらに、非正規化数というのもありますが、無視します。

*5:というのは嘘で、\(yw\) の符号に気をつけましょう。

*6:反例を見つけたら教えてください。

*7:binary64 は無限大たちを含むので \([-\infty, \infty]\) かも。

*8:古いのは確認していないのですが、当時はなかったらしい です? たぶん。

α(m, n) のお勉強 + 定数時間 decremental neighbor query

\(\alpha(n)\) や \(\alpha(m, n)\) というのがありますね。union-find の計算量解析で出てくる?らしい、ほぼ定数?みたいな関数?だと思われていがちなやつです。

一度お勉強したつもりですが、実はなにもわかっていなかったのでまたやりました(別にリンク先は読まなくてもいいです)。かなり長いので、暇なときにちょっとずつかいつまんで読むくらいにしてくれたらうれしいです。

rsk0315.hatenablog.com

わかっていない気がしていたこと:

  • \(\alpha(n)\) とはなにか + どんな性質を持つか
  • Ackermann 関数 \(A_m(n)\) とはなにか
    • \(\alpha(n)\) と \(A_m(n)\) の関係はどんなものか
  • \(\alpha(m, n)\) とはなにか + どんな性質を持つか
  • 計算量解析の文脈で、\(\alpha(n)\) や \(\alpha(m, n)\) はどういう過程で出てくるか

それから、これらの解析を用いて decremental neighbor query を \(\langle O(n), O(1)\rangle\) (amortized) で処理できるデータ構造を作れるので、それの話も書きます。何度か別で言及したりしていますが、このデータ構造によって線形時間でオフラインで LCA を処理できます。

これらを順を追って書いていきます。特に前半はやや天下りに感じますが許してください。

\(\alpha\) の性質に関する話

Ackermann 関数の逆関数として言及されがちですが、ここでは Ackermann 関数は介さずに定義していきます*1。Ackermann 関数との関係は後で示します。

Inverse Ackermann hierarchy

まずは \(\alpha(n)\) を定義する前に、inverse Ackermann hierarchy と呼ばれるものを定義します。 これは、\(\angled{\alpha_1(n), \alpha_2(n), \dots}\) という関数の(無限)列で、各関数の定義は以下です。

  • \(\alpha_1(n) = \ceil{n/2}\) for \(n\ge 1\),
  • \(\alpha_k(1) = 0\) for \(k\gt 1\),
  • \(\alpha_k(n) = 1 + \alpha_k(\alpha_{k-1}(n))\) for \(k\gt 1, n\gt 1\).
    • \(\alpha_k(n)\) は、\(\alpha_{k-1}\) を \(n\) に適用するのを繰り返して \(1\) にするための回数。
\(n\) 1 2 3 4 5 6 7 8 9 ... 15 16 17 ...
\(\alpha_1(n)\) 1 1 2 2 3 3 4 4 5 ... 8 8 9 ...
\(\alpha_2(n)\) 0 1 2 2 3 3 3 3 4 ... 4 4 5 ...
\(\alpha_3(n)\) 0 1 2 2 3 3 3 3 3 ... 3 3 4 ...

たとえば、\(\alpha_2(n) = \ceil{\log_2(n)}\) となります。また、\(\alpha_3(n)\) は反復対数 (iterated logarithm) と呼ばれる関数で、\(\log^{\star}(n)\) と書かれます。 より一般に、\(f\) を \(n\) に適用するのを繰り返して \(1\) にするための回数を表す関数を \(f^{\star}(n)\) と書くようなので、\(\alpha_k(n) = {\alpha_{k-1}}^{\star}(n)\) と書けそうです。

二変数バージョンであるところの \(\alpha(m, n)\) がもっと下で定義されるのですが、これとは別物なので注意してください。

さて、こうして定義される \(\alpha_k(n)\) が持つ性質を挙げていきます。

Claim 1: 任意の \(k\ge 1\) に対して以下が成り立つ。

  • \(\alpha_k(2) = 1\),
  • \(\alpha_k(3) = \alpha_k(4) = 2\),
  • \(\alpha_k(5) = \alpha_k(6) = 3\).

Proof: クリックして展開(以下同様)

\(k\) に関する帰納法で示す。

\(\alpha_1(2) = 1\) は明らか。 \(\alpha_k(2) = 1\) とすると、定義から \[ \begin{aligned} \alpha_{k+1}(2) &= 1 + \alpha_{k+1}(\alpha_k(2)) \\ &= 1 + \alpha_{k+1}(1) \\ &= 1 + 0 \\ &= 1. \quad\qed \end{aligned} \]

\(\alpha_k(3)\) から \(\alpha_k(6)\) についても同様に示せる。\(\qed\)

Claim 2: 任意の \(k\ge 1\) に対し、\(n\ge 4\implies \alpha_k(n) \le n-2\)。

Proof

\( (k, n)\) に関する帰納法で示す。 \(k = 1\) のときは明らかなので、以下 \(k\ge 2\) とする。\(4\le n\le 6\) については Claim 1 より成り立つ。

\( (k, n)\) より辞書順で小さいものについては成り立つとすると、 \[ \begin{aligned} \alpha_k(n) &= 1 + \alpha_k(\alpha_{k-1}(n)) \\ &\le 1 + \alpha_k(n-2) \\ &\le 1 + ( (n-2)-2) \\ &\lt n-2. \quad\qed \end{aligned} \]

Claim 3: 任意の \( (n, k)\) について \(\alpha_{k+1}(n)\le\alpha_k(n)\) が成り立つ。

さらに、\(k\ge 2\) のとき \(\alpha_k(n)\ge 4 \iff \alpha_{k+1}\lt \alpha_k(n)\)。

Proof

\(\alpha_k(n)\le 3\) のときは Claim 1 と同様にして \(\alpha_{k+1}(n)\le \alpha_k(n)\) が示せる。

\(\alpha_k(n)\ge 4\) のとき、Claim 2 から \[ \begin{aligned} \alpha_{k+1}(n) &= 1 + \alpha_{k+1}(\alpha_k(n)) \\ &= 1 + \alpha_k(n) - 2 \\ &\lt \alpha_k(n). \quad\qed \end{aligned} \]

Note: \(\alpha_2(1) = 0 \lt 1 = \alpha_1(1)\) なので、\(k=1\) では \(\Longleftarrow\) が成り立たなさそう。

Claim 4: \(k\ge 2 \implies \alpha_k(n) = o(n)\)*2

Proof

\(\alpha_2(n) = \ceil{\log_2(n)} = o(n)\) と \(\alpha_{k+1}(n)\le \alpha_k(n)\) から従う。\(\qed\)

Claim 5: \(k\ge 1 \implies \alpha_{k+1}(n) = o(\alpha_k(n))\)。

さらに、任意の正整数 \(r=O(1)\) に対して \(\alpha_{k+1}(n) = o({\alpha_k}^{(r)}(n))\)*3

Proof

Claim 4 より、 \[ \begin{aligned} \alpha_{k+1}(n) &= 1 + \alpha_{k+1}(\alpha_k(n)) \\ &= o(\alpha_k(n)). \quad\qed \end{aligned} \]

\(r\) についても \[ \begin{aligned} \alpha_{k+1}(n) &= r + \alpha_{k+1}({\alpha_k}^{(r)}(n)) \\ &= o({\alpha_k}^{(r)}(n)). \quad\qed \end{aligned} \]

\(\alpha(n)\) の話

本題ひとつめです*4

Claim 1, 3 より、任意の \(n\ge 5\) に対して \[ \alpha_1(n), \alpha_2(n), \alpha_3(n), \dots \] で定義される列は、最初は狭義単調減少し、その後 \(3, 3, 3, \dots\) となります。 たとえば、\(n = 9876!\) とすると、 \[ 3.06\times 10^{35163}, 116812, 6, 4, 3, 3, 3, \dots \] となります*5。 そこで、\(\alpha(n)\) を「何項目で \(3\) 以下になるか?」というのを返す関数として定義します。 \[ \alpha(n) \coloneqq \min\,\{k\mid \alpha_k(n)\le 3\}. \] 上記の例でいえば、\(\alpha(9876!) = 5\) となります。

いくつか値を示しておきます。

\(i\) \(\alpha(j) = i\) なる最小の \(j\)
1 1
2 7
3 9
4 17
5 65537

さて、こうして定義した \(\alpha(n)\) が持つ性質を挙げていきます。

Claim 6: 任意の定数 \(k\) に対して、\(\alpha(n) = o(\alpha_k(n))\)。

Proof

\(m = \alpha_{k+1}(n)\) とする。このとき、列 \[ \alpha_{k+1}(n), \alpha_{k+2}(n), \alpha_{k+3}(n), \dots \] を考えると、Claim 3 より \(\alpha_{k+m-2}(n) \le 3\) となる。 よって、 \[ \begin{aligned} \alpha(n) &\le k+m-2 \\ &= k+\alpha_{k+1}(n)-2 \\ &= o(\alpha_k(n)). \quad\qed \end{aligned} \]

Ackermann 関数との関係

さて、\(\alpha(n)\) を定義していくつか性質をわかったので、「union-find の計算量解析で出てくる?らしい、ほぼ定数?みたいな関数?」のようなふわっとした認識よりは成長した気がします。 それはそれとして、「Ackermann 関数の逆関数」などと呼ばれているのが妥当なのか?というのはよくわかっていません。

そこで、まず Ackermann 関数 \(A_m(n)\) を定義します。

  • \(A_0(n) = n + 1\) for \(n\ge 0\),
  • \(A_m(0) = A_{m-1}(1)\) for \(m\gt 0\),
  • \(A_m(n) = A_{m-1}(A_m(n-1))\) for \(m\gt 0, n\gt 0\).
    • \(A_m(n) = A_{m-1}^{(n+1)}(1)\) となります。

これらから、以下のことがわかります。

  • \(A_1(n) = n + 2\),
  • \(A_2(n) = 2n+3 = 2(n+3) - 3\),
  • \(A_3(n) = 2^{n+3}-3\).

+3 とか -3 とかがちょっと嫌な感じがしますが、諦めます。 さて、次のことが言えます。

Claim 7: \(m\ge 2 \implies \alpha_{m-1}(A_m(n)+3) = n+3\)。

Proof

\( (m, n)\) に関する帰納法で示す。

まず、\(m = 2\) のとき \[ \begin{aligned} \alpha_1(A_2(n)+3) &= \alpha_1( (2\cdot(n+3)-3) +3) \\ &= \alpha_1(2\cdot(n+3)) \\ &= n+3 \end{aligned} \] より成り立つ。

次に、固定した \(m\) 未満では成り立つとして、\(n = 0\) のとき \[ \begin{aligned} \alpha_{m-1}(A_m(0)+3) &= \alpha_{m-1}(A_{m-1}(1)+3) \\ &= 1 + \alpha_{m-1}(\alpha_{m-2}(A_{m-1}(1)+3)) \\ &= 1 + \alpha_{m-1}(1+3) \\ &= 1 + \alpha_{m-1}(4) \\ &= 3 \end{aligned} \] より成り立つ。Claim 1 より \(\alpha_{m-1}(4) = 2\) に注意せよ。

\( (m, n)\) より辞書順で小さいものについて成り立つとすると、 \[ \begin{aligned} \alpha_{m-1}(A_m(n)+3) &= \alpha_{m-1}(A_{m-1}(A_m(n-1))+3) \\ &= 1 + \alpha_{m-1}(\alpha_{m-2}(A_{m-1}(A_m(n-1))+3)) \\ &= 1 + \alpha_{m-1}(A_m(n-1)+3) \\ &= 1 + (n-1) + 3 \\ &= n + 3. \quad\qed \end{aligned} \]

Claim 8: \(n\ge 2\implies\alpha_{n-1}(n+3) = 3\)。

Proof

\(2\le n \le 3\) のとき、Claim 1 より \(\alpha_1(5) = \alpha_2(6) = 3\)。

\(n\gt 3\) とし、\(\alpha_1(n+3) = \ceil{(n+3)/2}\) と Claim 3 より、 \[ \begin{aligned} \alpha_{n-1}(n+3) &\le \max\,\{3, \alpha_1(n+3) - (n-3)\} \\ &= 3. \quad\qed \end{aligned} \]

Claim 9: \(\alpha(A_1(1)+3) = 1\) および \(n\ge 2\implies \alpha(A_n(n)+3) = n+1\)。

Proof

\(A_1(1) = 3\) より \(\alpha(A_1(1)+3) = 1\)。

以下 \(n\ge 2\) とする。 Claim 1, 7, 8 より、 \[ \begin{aligned} \alpha_{n-1}^{(2)}(A_n(n)+3) &= \alpha_{n-1}(n+3) = 3; \\ \alpha_{n-1}^{(3)}(A_n(n)+3) &= \alpha_{n-1}(3) = 2; \\ \alpha_{n-1}^{(4)}(A_n(n)+3) &= \alpha_{n-1}(2) = 1. \end{aligned} \] よって、定義から \(\alpha_n(A_n(n)+3) = 4\) とわかり、 \[ \begin{aligned} \alpha_{n+1}(A_n(n)+3) &= 1 + \alpha_{n+1}(\alpha_n(A_n(n)+3)) \\ &= 1 + \alpha_{n+1}(4); \\ &= 1 + 2 = 3; \\ \alpha(A_n(n)+3) &= n+1. \quad\qed \end{aligned} \]

こうして、(たかだか定数の差を厳密には示していないものの)\(\alpha(n)\) は \(f(n) = A_n(n)\) の逆関数くらいのものであることがわかりました。

\(A_m(\alpha_{m-1}(n))\) や \(A_n(\alpha_n(n))\) みたいなのについては割愛します*6

\(\alpha(m, n)\) の話

天下り感が強い気もしますが、次のように定義します。 \[ \alpha(m, n) \coloneqq \min\,\{k\mid \alpha_k(n)\le 3 + m/n\}. \]

Claim 10: \(\alpha(m, n) \le \alpha(n)\)。

Proof

\(\alpha_k(n)\le 3 \implies \alpha_k(n)\le 3+m/n\) より従う。\(\qed\)

Claim 11: \(\alpha(n\cdot\alpha_i(n), n) \le i\)。特に、\(\alpha(n\ceil{\log_2(n)}, n) \le 2\)。

Proof

\[ \alpha_i(n)\le 3 + (n\cdot\alpha_i(n))/n) = 3 + \alpha_i(n) \] が自明に成り立つことから従う。後者は、\(i = 2\) の場合。\(\qed\)

注意書き

\(\alpha(n)\) は、\(A\) を使わずに定義する場合と \(A\) を使って定義する場合があります。 \(A\) を使う場合は \(\alpha(n) = \min\,\{k\mid A_k(1)\ge n\}\) のような感じの形で定義したりしますが、中身の右辺が \(n\) ではなく \(\log_2(n)\) とかだったり亜種があります。あるいは、\(A\) の定義自体にも亜種があります。 また、inverse Ackermann hierarchy で定義する場合も、\(\alpha_1(n)\) が \(\floor{n/2}\) だったり \(\ceil{n/2}\) だったり亜種があります。 各定義で高々定数の差はあったりはします*7

文脈に合わせて都合よく定義して、オーダーは同じなので ok!みたいな感覚でやっちゃっていい気がします*8

\(\alpha\) の計算量解析での話

ある程度 \(\alpha_k(n)\) や \(\alpha(n)\) と仲よしになったと思うので、ここからは計算量解析の文脈でどう出てくるかというのを書いていきます。 \(A_m(n)\) はもう御役御免です、ごめんね。

この辺は、参考文献に載せたスライドを読むとわかりやすいかもしれません。 数式の部分だけ読みたい人は数式の部分だけ探して読んでもいいかもしれません。

\(\alpha(n)\) が出てくるデータ構造

モノイドの区間積を答えるデータ構造を考えます*9。更新はないものとしてよいです。

まず、「前処理をなるべく少なく行いつつ、各クエリでは 1 回のモノイド演算で済むようにしたい」というのを考えます。

素数が 1 のときは自明なので、2 以上とします。 中間で半分に分けて、前半は「その要素から中間までの累積モノイド積」、後半は「中間からその要素までの累積モノイド積」を求めておきます。 これにより、中間をまたぐ区間のクエリについてはそれを用いて 1 回の演算で答えられます。 中間をまたがない区間のために、再帰的に各々処理します*10

各段でのモノイド演算の回数はたかだか \(n\) で、段数は \( (n/2)^{\star} = \log_2(n)\) なので、前処理は \(n\log_2(n)\) 回程度です。

「長さ \(n\) の配列のモノイド区間積に対し、各クエリを \(i\) 回のモノイド演算で済ませるための前処理の最小回数」を \(S_i(n)\) と書くことにすると、\(S_1(n) \le n\log_2(n)\) と書けます。

次に、このデータ構造を用いて、クエリあたり 3 回で済ませるデータ構造を考えます。

  • 配列を \(\log_2(n)\) 個ずつのブロック \(n/\log_2(n)\) 個に分ける
    • 各ブロックについて、前から・後ろからの累積モノイド積を求めておく
    • 高々 \(2n\) 回
  • ブロック全体のモノイド積を 1 要素と見たときの、区間クエリの前処理をする
    • 配列サイズが \(n/\log_2(n)\) となる
    • 高々 \(S_1(n/\log_2(n)) \le n\) 回
  • ブロック内のクエリを 3 回で済ませられるようにする。
    • 再帰的に行い、\(n/\log(n)\cdot S(\log_2(n))\) 回

再帰の段数は \(\log_2^{\star}(n)\) となり、\(S_3(n)\le 3n\log_2^{\star}(n)\) となります。

同様の手続きを繰り返し、\(S_{2k-1}(n)\le (2k-1)\cdot n\cdot f(n)\) を達成するデータ構造を用いて \(S_{2k+1}(n)\le (2k+1)\cdot n\cdot f^{\star}(n)\) を達成できます。すなわち、以下のようになります。 \[ S_{2k+1}(n) = (2k+1)\cdot n\cdot\log\overbrace{{}_{\!2}^{\star\star\,\dots\,\star}}^{k\text{ times}}(n). \] この \(\log^{\star\star\,\dots\,\star}\) の部分を定数にするには \(k\) をいくらにすればよいでしょう。というのを考えると、次のような関数が欲しくなります。 \[ \alpha(n) = \min\,\{k\mid \log\overbrace{{}_{\!2}^{\star\star\,\dots\,\star}}^{k\text{ times}}(n)\le 2\}. \] 前述の \(\alpha\) と違って \(\ceil{\;}\) をしていなかったりする関係で定数部分が違ったりしますが、同じようなものです。こうした過程で出てくるわけですね。

すなわち、\(S_{2\alpha(n)+1} \le (2\cdot\alpha(n)+1)\cdot n = O(n\cdot\alpha(n))\) となります*11

\(\alpha(m, n)\) が出てくるデータ構造

みんな大好き union-find。

ちょっと疲れちゃったので、導出は抜きにして数式だけ書きます。 図や詳細などが必要な人は、同じくリンク先のスライドを見てください。

「頂点数 \(n\) かつランクの最大値 \(r\) の rank forest に対して操作を任意に \(m\) 回行うときの最大コスト」を \(f(m, n, r)\) と定義します。

さて、 \[ f(m, n, r) \le k\cdot m+2n\cdot g(r) \implies f(m, n, r) \le (k+1)\cdot m+2n\cdot g^{\diamond}(r) % ({\ceil{\log_2}}\circ g)^{\star}(r) \] が言えるらしいです。ここで、\(g^{\diamond} = ({\ceil{\log_2}}\circ g)^{\star}\) です。あるいは次のように書けます。 \[ g^{\diamond}(r) = \begin{cases} 0, & \text{if }r\le 1; \\ 1 + g^{\diamond}(\ceil{\log_2(g(r))}), & \text{otherwise}. \end{cases} \] \(m\) 側の項をちょっと悪くして \(n\) 側の項をめちゃよくするみたいなイメージだと思います。

自明な上界として \(f(m, n, r) \le (r-1)\cdot n\) があるらしい*12ことから、\(g_0(r) = \ceil{(r-1)/2}\) とします*13。また \(r\le\floor{\log_2(n)}\) を踏まえて「\(n\) 側の項を定数にするには?」というのを考えると、次のようになります。 \[ \alpha(n) = \min\,\{k\mid g\,\overbrace{{}_{\!0}^{\diamond\diamond\,\dots\,\diamond}}^{k\text{ times}}(\floor{\log_2(n)})\le 2\}. \] あるいは、「\(n\) 側の項を \(O(m)\) にするには?」というのを考えると、次のようになります*14。 \[ \alpha(m, n) = \min\,\{k\mid g\,\overbrace{{}_{\!0}^{\diamond\diamond\,\dots\,\diamond}}^{k\text{ times}}(\floor{\log_2(n)})\le 2+m/n\}. \] こうすれば、 \[ \begin{aligned} f(m, n, r) &\le \alpha(m, n)\cdot m+2n\cdot (2+m/n) \\ &\le \alpha(m, n)\cdot m + 4n+2m \end{aligned} \] となり、操作が \(m\) 回なので、操作 1 回あたり \(O(\alpha(m, n))\) となります*15

こちらも前述の \(\alpha(m, n)\) とは若干違いますが、大した影響ではないでしょう。\(\alpha_k = {\alpha_{k-1}}^{\star}\) で定義したものと \(\alpha_k = {\alpha_{k-1}}^{\diamond}\) で定義したものでオーダーが変わるかは考えていません。\(\ceil{\log{}\circ}\) をつけた程度では \(\alpha\) に影響しないのではないかという気がしています。Claim 5 あたりを使うとよいかもしれません。

Decremental neighbor query

さて、後半です。えびちゃんは疲れてきました。記事を分けるべきだったかもしれません。

問題設定

  • \(A \gets \{0, 1, \dots, n-1\}\) で初期化。
  • 与えられた \(i\) に対して \(A \gets A \setminus \{i\}\) で更新。
  • 与えられた \(i\) に対して \(\min\,\{j\in A\mid j\ge i\}\) を出力。
  • 与えられた \(i\) に対して \(\max\,\{j\in A\mid j\le i\}\) を出力。
    • Note: \(i\in A\) とは限らない。

後者二つ (successors query, predecessor query) を合わせて neighbor query と呼びます。 要素が減るだけなので decremental です。初期化は任意の \(S\subseteq \{0, 1, \dots, n-1\}\) を与えてもよいです。

解法

熨斗袋先生ありがとう*16

word size \(w \ge \log_2(n)\) で考えます。 長さ \(n/\log_2(n)\) の配列を取り、各要素は \(\log_2(n)\) bits です。 neighbor query の答えが、\(i\) を含む要素と同じ要素にあるときは、\(O(1)\)-time MSB やビット演算などを用いて定数時間でできます*17

そうでないときのこと考えます。隣り合う要素に対して、それらの要素の bit がすべて 0 になったとき、union-find でそれらをつなぎます。 根を管理するのと同様にして、「その連結成分(区間)の最左はどこか?」「最右はどこか?」というのも管理しておきます。 その位置を union-find で求めた後は、上記同様に要素内で neighbor query が完結するので、定数時間で答えられます。

さて、union-find の計算量について考えます。 この union-find の要素数は \(n/\log(n)\) なので、\(n \ge n/\log(n)\cdot \log(n/\log(n))\) 回クエリをすると一回あたり \(\alpha(n, n/\log(n)) \le 2 = O(1)\) 時間です。 \(n\) 回操作をしたとき \(O(n)\) 時間になるので、\(n\) 回操作をする前のコストに関しては初期化のコストにでも押しつけておけばよいです。

これにより、\(\angled{O(n), O(1)}\) (amortized) で decremental neighbor query が解けたことになります。

これを使って解ける問題たち:

参考文献

そういえば、元ツイ時点でのえびちゃんは、配列で(dancing links をやるときみたいに)連結リストを持っていたので、neighbor query には答えられませんでした。

こっちは論文たちです。

  • Tarjan, Robert Endre. "Efficiency of a good but not linear set union algorithm." Journal of the ACM (JACM) 22, no. 2 (1975): 215–225.
    • union-find の計算量解析の元論文
  • Gabow, Harold N., Zvi Galil, Thomas Spencer, and Robert E. Tarjan. "Efficient algorithms for finding minimum spanning trees in undirected and directed graphs." Combinatorica 6, no. 2 (1986): 109–122.
    • \(\beta(m, n) = \min\,\{i\mid\log^{(i)}(n)\le m/n\}\) で定義される関数が出てくる

ツイート振り返りのコーナー

こんなツイートがきっかけでこの記事を書くことになったらしいです。信じられません。

そういえば、\(\arcsin\) を先に定義して、その逆関数として \(\sin\) を定義する話もあったような気がします。

これは別件ですが、これはこれで面白いです。

帰納法で詰まっていたときのツイートです。

おまけ

ここでは記事中で最初の \(\alpha_k(n)\) や \(\alpha(n)\) に関する定義を採用するものとして、Claim たちをいくつか書きます。あんまり詰めてません。

Claim 12: \(\alpha_k(n+1) - \alpha_k(n) \le 1\)。

Proof

\( (k, n)\) の帰納法で示す。

\(k = 1\) および \(\alpha_k(n)\le 3\) については成り立つ。 \( (k, n)\) 未満について成り立つとすると、 \[ \begin{aligned} \alpha_k(n+1) - \alpha_k(n) &= (1 + \alpha_k(\alpha_{k-1}(n+1))) - (1 + \alpha_k(\alpha_{k-1}(n))) \\ &= \alpha_k(\alpha_{k-1}(n+1)) - \alpha_k(\alpha_{k-1}(n)) \\ &\le \alpha_k(\alpha_{k-1}(n)+1) - \alpha_k(\alpha_{k-1}(n)). \end{aligned} \] \(\alpha_{k-1}(n) = 3\) であれば \(\alpha_k(4)-\alpha_k(3)=0\le 1\) より成り立つ。 そうでないとき、Claim 3 より \(n\lt\alpha_{k-1}(n)\) なので、帰納法の仮定より成り立つ。\(\qed\)

\(\alpha_k(n+\bullet)-\alpha_k(n)\le 1\) の \(\bullet\) の部分をもっと大きくできる? 端数は適当だけど \(A_k(2)-A_k(1)\) かなんかそんな感じに。

Claim 13: 非負整数 \(r\) に対して \(\alpha_{k+1}(n+r)-\alpha_{k+1}(n) \le \alpha_k(n+r)-\alpha_k(n)\)。

Proof

\[ \alpha_{k+1}(n+r)-\alpha_k(n+r) \le \alpha_{k+1}(n)-\alpha_k(n) \] より、\(f_k(n) = \alpha_{k+1}(n) - \alpha_k(n)\) が \(n\) に対して単調減少であることを示す。 \[ \begin{aligned} f_k(n) &= \alpha_{k+1}(n) - \alpha_k(n) \\ &= 1 + \alpha_{k+1}(\alpha_k(n)) - \alpha_k(n). \end{aligned} \] \(m = \alpha_k(n)\) とすると、Claim 12(と \(m\) の単調性)から \[ \begin{aligned} 1 + \alpha_{k+1}(m + 1) - (m + 1) &\le 1 + 1 + \alpha_{k+1}(m) - (m + 1) \\ &= 1 + \alpha_{k+1}(m) - m \end{aligned} \] より従う。\(\qed\)

Claim 14: \(k\ge 2 \implies \alpha_k(2n)-\alpha_k(n)\le 1\)。

Proof

\(\ceil{\log_2(2n)}-\ceil{\log_2(n)} = 1\) と Claim 13 より従う。

Claim 15: \(\alpha(\alpha_k(n)) = \Theta(\alpha(n))\)。

Proof

\(n\) は十分大きいとする。

\(\alpha_{k+1}(n) - \alpha_{k+1}(\alpha_k(n)) = 1\) より、Claim 12 を繰り返し適用したりして \(\alpha_{\alpha(n)-1}(\alpha_k(n)) \le 4\) になって、\(\alpha(\alpha_{k+1}(n)) \in\{\alpha(n)-1, \alpha(n)\}\) とかになりそう。\(\qed\)

Claim 16: \(m\ge 3, n\ge 1 \implies \alpha_{m-2}(A_m(n)+3) = A_m(n-1)+3\) および \(m\ge 3 \implies \alpha_{m-2}(A_m(0)+3) = 4\)。

Proof Claim 7 と \(A_m\) の定義から \[ \begin{aligned} \alpha_{m-2}(A_m(n)+3) &= \alpha_{m-2}(A_{m-1}(A_m(n-1))+3) \\ &= A_m(n-1)+3. \end{aligned} \] また、 \[ \begin{aligned} \alpha_{m-2}(A_m(0)+3) &= \alpha_{m-2}(A_{m-1}(1)+3) \\ &= 1+3 = 4.\quad\qed \end{aligned} \]

Claim 17: \(m\ge 1, n\ge 0 \implies A_m(n) \ge n + 2\)。

Proof

\( (m, n)\) に関する帰納法で示す。

\(m = 1\) のとき定義より等号が成り立つ。

固定した \(m\) 未満で成り立つとすると、\(A_m(0) = A_{m-1}(1) \ge 1 + 2 \ge 0 + 2\) より成り立つ。

\( (m, n)\) より辞書順で小さいものについて成り立つとすると、 \[ \begin{aligned} A_m(n) &= A_{m-1}(A_m(n-1)) \\ &\ge A_m(n-1) + 2 \\ &\ge n+1 + 2 \\ &\ge n+2. \quad\qed \end{aligned} \]

Claim 2 と似ている。

Claim 18: \(m\ge 2, n\ge 1 \implies A_m(n)-2 \ge A_m(n-1)\)。

Proof

帰納法で示す。

\(m = 2\) のとき \(A_2(n)-2 = 2n+1 = A_2(n-1)\) より成り立つ。

固定した \(m\) 未満で成り立つとすると、Claim 17 より \[ \begin{aligned} A_m(n) - 2 &= A_{m-1}(A_m(n-1)) - 2 \\ &\ge (A_m(n-1)+2) - 2 \\ &= A_m(n-1). \quad\qed \end{aligned} \]

Claim 19: 固定した正整数 \(i = O(1)\) に対して、\(\alpha'_i = \alpha_i\) および \(\alpha'_{i+1} = (\alpha_i\circ\alpha'_i)^{\star}\) とし、\(\alpha'(n) = \min\,\{k\mid \alpha'_k(n)\le 3\}\) とする。このとき、\(\alpha'(n) = \Theta(\alpha(n))\)。

Proof

Note: \(i = 2\) のとき、上での \(g_0^{\diamond}\) に対応する。

\(\alpha'_m(n) \le \alpha_m(n)\) は明らか。

まず、\(\alpha'_{i+1} = (\alpha_i\circ\alpha_i)^{\star}\) なので、\(\alpha'_{i+1}(n) = \ceil{\alpha_{i+1}(n)/2}\)。

\[ \begin{aligned} \alpha_{i+2}(A_{i+3}(n)+3) &= 1 + \alpha_{i+2}(\alpha_{i+1}(A_{i+3}(n)+3)) \\ &= 1 + \alpha_{i+2}(A_{i+3}(n-1)+3) \\ &= n + \alpha_{i+2}(A_{i+3}(0)+3) = n+3 \end{aligned} \] であるが、 \[ \begin{aligned} \alpha'_{i+2}(A_{i+3}(n)+3) &= 1 + \alpha'_{i+2}( (\alpha_i\circ\alpha'_{i+1})\,(A_{i+3}(n)+3)) \\ &= 1 + \alpha'_{i+2}(\alpha_i(\alpha'_{i+1}(A_{i+3}(n)+3))) \\ &= 1 + \alpha'_{i+2}(\alpha_i(\ceil{\alpha_{i+1}(A_{i+3}(n)+3)/2})) \\ &\ge 1 + \alpha'_{i+2}(\alpha_i(2\ceil{\alpha_{i+1}(A_{i+3}(n)+3)/2})-1) \\ &\ge 1 + \alpha'_{i+2}(\alpha_i(\alpha_{i+1}(A_{i+3}(n)+3))-1) \\ &\ge 1 + \alpha'_{i+2}(\alpha_i(A_{i+3}(n-1)+3)-1) \\ &= 1 + \alpha'_{i+2}(\alpha_i(A_{i+2}(A_{i+3}(n-2))+3)-1) \\ &= 1 + \alpha'_{i+2}(A_{i+2}(A_{i+3}(n-2)-1)+3-1) \\ &\ge 1 + \alpha'_{i+2}(A_{i+2}(A_{i+3}(n-2)-2)+3) \\ &\ge 1 + \alpha'_{i+2}(A_{i+2}(A_{i+3}(n-3))+3) \\ &= 1 + \alpha'_{i+2}(A_{i+3}(n-2)+3) \\ &\ge \begin{cases} \frac{n}{2} + \alpha'_{i+2}(A_{i+3}(0)+3), & \text{if }n\equiv 0\pmod{2}; \\ \frac{n-1}{2} + \alpha'_{i+2}(A_{i+3}(1)+3), & \text{otherwise}. \end{cases} \end{aligned} \] ここで、 \[ \begin{aligned} \alpha'_{i+2}(A_{i+3}(0)+3) &\ge 1; \\ \alpha'_{i+2}(A_{i+3}(1)+3) &= 1 + \alpha'_{i+2}(\alpha_i(\ceil{\alpha_{i+1}(A_{i+3}(1)+3)/2})) \\ &= 1 + \alpha'_{i+2}(\alpha_i(\ceil{(A_{i+3}(0)+3)/2})) \\ &\ge 1 + \alpha'_{i+2}(\alpha_i(2\ceil{(A_{i+3}(0)+3)/2})-1) \\ &\ge 1 + \alpha'_{i+2}(\alpha_i(A_{i+3}(0)+3)-1) \\ &= 1 + \alpha'_{i+2}(\alpha_i(A_{i+2}(1)+3)-1) \\ &= 1 + \alpha'_{i+2}(A_{i+2}(0)+3-1) \\ &= 2 \end{aligned} \] なので、 \[ \begin{aligned} \alpha'_{i+2}(A_{i+3}(n)+3) &\ge \Floor{\frac{n+3}{2}} \\ &= \Floor{\frac{\alpha_{i+2}(A_{i+3}(n)+3)}{2}} \end{aligned} \] が言える。

同様にして、\(m\gt i+2\) についても \(\alpha'_m(A_{m+1}(n)+3) \ge \floor{\alpha_m(A_{m+1}(n))/2}\) が示せる。 これより、\(\alpha'(n) = \Theta(\alpha(n))\) も従う。\(\qed\)

もっとちゃんと見積もればもっと綺麗に示せそう。

Claim 20: \(k\ge 2\) に関して、\(\log_2(n)\) に関するいい感じの不等式が成り立つならば、その不等式の \(\log_2\) を \(\alpha_k\) に置き換えた不等式も成り立つ。

これは Claim か?

おきもち

Claim 14 の「\(\log_2\) で成り立つので Claim 13 から従う」の一般化みたいなことがしたかったです。Claim 14 の引数部分を \(A_m\) とか \(\alpha_m\) とかにしたいのかも。

このレベルの \(\alpha_k\) でもこの処理を無効化できるのだからもっと上のレベルなら当然無効でしょみたいな。

おわり

にゃんこです。つかれましたね。

*1:サンタさんとは独立にサンタさんの帽子を定義できるという例え話があります。

*2:ざっくり言って、\(o(f(n) )\) はオーダーが \(f(n)\) より真に小さいことを表す。

*3:\(f^{(r)}(n)\) は \(f\) を \(n\) に \(r\) 回適用したものを表す。\(f^{(0)}(r) = r\) と \(f^{(r)}(n) = f(f^{(r-1)}(n) )\)。\(r\) 階微分ではないです。

*4:残りの本題がどこだと思うかは読む人の主観に任せます。

*5:\(9876!\) がめちゃ大きいのに対し、すぐ \(3, 3, \dots\) が現れるなぁというのが実感できます。

*6:\(n\) が \(\alpha\) によって小さく丸められた後 \(A\) によって大きく飛ばされるので、\(n\) 自体とはそれなりに差が出てきそう。上からよしなに抑えられるとかの関係を示すのかなぁ。

*7:これを読んでいる人も定数倍を気にしない人が多いでしょうから、差なんてさらに気にしないことと思います。

*8:定数の差しかないことはちゃんと示した方がいいと思います。

*9:RMQ みたいなのだと思ってください。任意の演算で。

*10:disjoint sparse table ですね。

*11:「\(\alpha(n)\) は定数」学派の人的には、任意のモノイド積の区間クエリを、線形時間前処理・定数時間クエリ処理できるわけですね。まぁ「log は定数」学派の人はセグ木でそれができるんですけど。ちなみに、ちょっと工夫するだけで \(\angled{O(n), O(\alpha(n) )}\) にできます。

*12:trivial bound と書かれていますが、今ぼーっとしていてよくわかりません。疲れてきました。

*13:\(k=0\) として \(g(r) = g_0(r)\) とすればいい感じになるはず。

*14:スライドでは \(1+m/n\) になっていますが、\(1\) で問題ないかはあまりわかっていません。どうせ定数の差だと思うのですが。あるいは \(m\ge n\) を仮定しているかも。

*15:\(n\) に関する項は、初期化のコストにでも押しつければいいと思います。

*16:えびちゃんのツイートに対して熨斗袋先生が言ってくれたデータ構造なので思い入れが深いです。

*17:あるいは、leading zero を求める演算などがマシン語命令にあるので \(O(1)\) 時間でできるという仮定にしてもよいと思います。

遅延セグ木を思いつく流れを追う用の記事

敬遠されがちみたいな印象がある[要出典]*1ので、書いてみます。

問題

以下のような問題設定を考えます。

集合 \(S\), \(T\) と、それに関する演算 \(\circ: S\times S\to S\), \(\star: S\times T\to S\) があります。 変数 \(a_1, a_2\in S\) があり、初期値が与えられます。 あなたは、変数 \(a'_1, a'_2, a' \in S\) と \(x' \in T\) のみ読み書きできます。 以下のクエリを処理してください。

  • 単一クエリ
    • \(x\) が与えられ、\(a_1 \gets a_1 \star x\) で更新する。
    • \(x\) が与えられ、\(a_2 \gets a_2 \star x\) で更新する。
    • \(a_1\) を答える。
    • \(a_2\) を答える。
  • 複数クエリ
    • \(x\) が与えられ、\(a_1 \gets a_1 \star x\) および \(a_2 \gets a_2 \star x\) で更新する。
    • \(a_1 \circ a_2\) を答える。

ただし、複数クエリが与えられたときは、\(a'\) と \(x'\) にのみアクセスできるとします。

単一クエリにおいては、アクセスは自由とします。 たとえば、\(a_1\) を答えるときに \(a'_2\) や \(x'\) を更新してもよいとします。

ここでは \(\circ\), \(\star\) が持つ性質には特に触れていません。 これ以降で問題を解きながら「こういう性質があれば解けるな〜」というのを探していく流れになります。

記号が苦手な人は、\(S\), \(T\) が int、\(\circ\) が min、\(\star\) が + のような例をイメージするといいかもしれません*2

観察 1

まず、\(a_1 \circ a_2\) を答える際に \(a'_1\), \(a'_2\) に触れられないので、\(a' = a_1 \circ a_2\) にしておきたいな〜という気持ちになります。

それから、更新する際にもそれらに触れられないので、\( (a_1 \star x) \circ (a_2 \star x) = (a_1\circ a_2)\star x\) みたいな性質が成り立ってほしいな〜となります。前者と合わせて \(a \gets a\circ x\) と更新してよくなるためです。

さらに、このとき、実際には \(a_1 = a'_1 \star x\) のように乖離が生じるため、「\(a'_1\gets a'_1 \star x\) する必要がある」というのを覚えておく必要があります。これを \(x'\) で管理しておきます。

また、連続して複数更新クエリ(値はそれぞれ \(x\), \(y\) とする)が来ると、\(T\) に関する変数は一つしか持っていないので困ります。 そこで、ある演算 \(\ast: T\times T\to T\) があり、\(x'\gets x\ast y\) としてもよいという性質があるとうれしいです。 すなわち、一回ずつ更新する \( (a'_1\star x)\star y\) と、まとめて更新する \(a'_1 \star (x\ast y)\) が同じになってほしいです。 \(a'_2\) についても同様で、一般に \(S\) の要素全部についてそうであってほしいです。

欲しい性質まとめ

ここまでで欲しくなった性質は次の通りです。

  • \( (a_1 \star x)\circ(a_2 \star x) = (a_1\circ a_2) \star x\)
  • \( (a \star x)\star y = a \star (x\ast y)\) を満たす \(\ast: T\times T\to T\) が存在する

観察 2

上記の性質が成り立つとして、方針を詰めていきます。

以下の不変条件が成り立つように管理します。

  • \(a_1 = a'_1 \star x'\)
  • \(a_2 = a'_2 \star x'\)
  • \(a' = a_1 \circ a_2\)

これらより、\(a' = (a'_1 \circ a'_2) \star x'\) が成り立っていることになります。

複数クエリに関しては、\(a'\) を答えるのと、\(a' \gets a'\star x\) および \(x' \gets x'\ast x\) で更新するだけなので、単一クエリについて考えます。

単一クエリでは、変数に自由にアクセスできるので、\(a'_1 = a_1\) が成り立つようにします。\(a'_1 \gets a'_1 \star x'\) です。\(a'_2\) についても同様です。 この際、不変条件から、\(a_1 = a_1 \star x'\) になってほしいので、そういう元 \(\mathrm{id}_T\) があって \(a_1 = a_1\star \mathrm{id}_T\) が成り立ってほしいです*3。これを使って \(x'\gets \mathrm{id}_T\) とします。

もう \(a'_1 = a_1\) になったので、出力クエリでは \(a'_1\) を答えればよいです。更新クエリにおいては、\(a'_1 \gets a'_1 \star x\) で更新し、\(a'\gets a'_1\circ a'_2\) と直しておきます。こちらも \(a'_2\) についても同様です。

\(a'_1 = a_1\) にすることを、\(a'_1\) の遅延を解消する (force) と呼ぶことにします*4

本題

ここまで来れば、遅延セグ木は実質わかったと言えそうです。

前問で 2 つだけの要素 \(a_1, a_2\in S\) だったのを、4 つの要素 \(a_1, a_2, a_3, a_4\in S\) にしてみます。

           [ (a[1234], x[1234]) ]
 [ (a[12], x[12]) ]      [ (a[34], x[34]) ]
[ a[1] ]    [ a[2] ]    [ a[3] ]    [ a[4] ]

アクセスできる変数は、関連するノードの(真の (proper))祖先の直接の子のみとします。たとえば、a[3] に関するクエリが来た場合は a[3] a[4] a[12] x[12] a[34] x[34] a[1234] x[1234]a[12] に関するクエリが来た場合は a[12] x[12] a[34] x[34] a[1234] x[1234]a[1234] に関するクエリが来た場合は a[1234] x[1234] です。

こんな感じでイメージすると、前問の

    [ (a, x) ]         [ (a[12], x[12]) ]
[ a[1] ]  [ a[2] ] =~ [ a[1] ]    [ a[2] ]

    [ (a, x) ]         [ (a[34], x[34]) ]
[ a[1] ]  [ a[2] ] =~ [ a[3] ]    [ a[4] ]

    [ (a, x) ]                [ (a[1234], x[1234]) ]
[ a[1] ]  [ a[2] ] =~ [ (a[12], x[12]) ]  [ (a[34], x[34]) ]

との対応づけがわかりそうです。最後(木で言うと根)に関してのみ、子が \(x\) を持っているのが違いますが、\(x\) はそのまま子に渡してよい (\( (a_1\circ a_2)\star x = (a_1\star x)\circ(a_2\star x)\) から従う) ので、\(a'_{\text{child}}\) の更新の際に、加えて \(x'_{\text{child}} \gets x'_{\text{child}} \ast x'_{\text{parent}}\) で更新すればよいです。

\(a_1\circ\dots\circ a_4\) を求めるクエリにおいては、(適宜 force して更新してから)ふつうのセグ木のようにまとめていくので、\( (S, \circ)\) はモノイドになっていてほしいです。

また、\(T\) に関しても、\(x_1 \ast x_2\) で更新されてから \(x_3\) で更新された場合と、\(x_1\) で更新されてから \(x_2\ast x_3\) で更新された場合とで値が変わると困るので、\( (T, \ast)\) もモノイドになってほしいです。

これらより、必要なノードにアクセスするだけでクエリ処理できることになりました。 もちろん、4 つではなくて一般に \(n\) 個について適用でき、必要なノードは \(O(\log(n))\) 個なので、\(O(\log(n))\) 時間で処理できます。

成り立ってほしい性質

  • \( (S, \circ)\) はモノイドを成す*5
  • \( (a_1 \star x)\circ(a_2 \star x) = (a_1\circ a_2) \star x\)
  • \( (a \star x)\star y = a \star (x\ast y)\) を満たす \(\ast: T\times T\to T\) が存在する
  • \( (T, \ast)\) はモノイドを成し、単位元が \(\mathrm{id}_T\)
  • \(a \star \mathrm{id}_T = a\)

その他

遅延セグ木で問題を解く際には、\(\ast\) を高速に計算できるように考えたり、\(\star\) に関する性質が成り立つように \(S\) や \(T\) を工夫することが基本になりそうです。

初心者が勘違いしがちですが、区間 \(\star\) 更新クエリと区間 \(\circ\)-fold クエリの内容によって、\(S\) が持つべき情報(あるいは定義の可否)が変わります。区間加算クエリと区間 \(\min\) クエリが可能だからといって、区間加算クエリと区間任意-fold クエリができたり、区間任意更新クエリと区間 \(\min\) クエリができたりするとは言えません。

また、区間加算クエリを区間和クエリと一緒にやるためには、\(S\) に要素数も持たせる必要もあります*6int とかのペアにするみたいな感じですね。

参考文献

qiita.com

opt-cp.com

ACL のドキュメント + Appendix も読むとよいです。

おわりねこ

遅延セグ木とか全方位木 DP とかのように、うまいことやっておいて汎用的に使えるライブラリを作る・考えるのは、よい勉強になると思います。遅延セグ木は ACL にもありますが、実装の考え方をわかっていると、どのような代数的構造が扱えるのかを丸暗記せずに済むのでよさそうです。

抽象化ライブラリの作成はコンテストで問題を解くよりも得てして大変なことが多いです。

競技プログラミングの勉強法としては不向きかもしれませんが、アルゴリズムプログラミング言語に対する理解を深める方法の 1 つとしてはありだと個人的に思っています。

atcoder.jp

競プロで出てくる抽象化ライブラリ、大抵がモノイド (+ なんか) が乗ってがち*7なので、競プロ er モノイドとは仲いいがち。BIT はモノイドではないけど。

ダブリングとかもモノイドが乗ってそう。木や functional graph 上でやるやつとか。

WM とか CHT とか slope trick とかは抽象化したりしにくそう。\(\langle O(n), O(1)\rangle\) decremental neighbor query とかも整数モノだからつらそう。

WM + なんかで重みつき矩形和みたいなのは可換モノイドが乗る? Union-find にポテンシャルをつけるやつ(最近あんまり見ない)には群が乗る?

文字列系のやつを除けば、初等的な競プロで出てくるデータ構造ってこのくらいしかない気がするね。

rsk0315.hatenablog.com

CHT も知りたい人は読んでね。

*1:別に個人の感想なんだから出典とかいらなくない?

*2:実際には int は + について閉じていないですが、そういうのは一旦無視です。

*3:\(\star\) しても \(S\) の値を変えないような元。\(\mathrm{id}_T\) は、\( (T, \ast)\) についても単位元になりそうです。実際にはなるとは限らない場合もあるかも。そんなことないかも。

*4:関数型言語とかの遅延評価の文脈でそういう用語が使われるみたいなので。遅延をやめて、評価することを強制するみたいな意味合いかも。

*5:モノイドをなすがすきなす

*6:素数に関しては、こうした考察から必要になってくるものなので、モノイド側が計算するのが自然で、セグ木側で渡してあげるのは不自然かなという気持ちがあります。区間代入クエリと区間ロリハクエリとかでは、要素数ではなく b素数 が必要になったりするので。

*7:foldable queue (SWAG)、tree catamorphism (全方位木 DP)、とか。

クエリ先読みなし任意順序 CHT

あまり知られていないような気がするので書きます。欲しくなる状況自体が珍しいという話かも。

適当に調べると、傾きの単調性を使って deque で処理するやつか、先読みして Li-Chao tree で処理するやつがたくさん見つかる気がします。 先読みなし・単調性の仮定なしで、\(O(\log(n))\) time per query でやります。

メモ

ここあんまり考えてません。

よくある log 段積んで静的データ構造を動的にするテクとかでもどうにかなる? 定数時間でクエリ処理するの自体つらそうだから log ふたつはついちゃったりしそう。

メモおわり。

定義

まず、解きたい問題を書いておきます。 以下のクエリを処理したいです。

  • \(S \gets \emptyset\) で初期化。
  • 与えられた \( (a, b)\) に対して \(S \gets S \cup \{(a, b)\}\) で更新。
  • 与えられた \(x\) に対して \(\min_{(a, b)\in S} ax+b\) を出力。

直感的には、「xy-平面に直線を追加する」「今ある直線のうち、特定の x 座標における y 座標の最小値を答える」というものです。

この形式のクエリによって高速化できる DP があるのですが、問題例のリンクだけいくつか貼って詳細は割愛します。 AtCoder, Codeforces, AOJ

CHT というのは convex hull trick の略です。どこ発祥の呼び方なんだろう。

アイディア

たとえば、\(y = 2x+3\), \(y = x+1\), \(y = 3\), \(y = -x+4\), \(y = -x+5\) の直線がある状況を考えてみます。

   x | -5 | -4 | -3 | -2 | -1 |  0 |  1 |  2 |  3 |  4 |  5 |
2x+3 | -7 | -5 | -3 | -1 |  1 |  3 |  5 |  7 |  9 | 11 | 13 |
 x+1 | -4 | -3 | -2 | -1 |  0 |  1 |  2 |  3 |  4 |  5 |  6 |
   3 |  3 |  3 |  3 |  3 |  3 |  3 |  3 |  3 |  3 |  3 |  3 |
-x+4 |  9 |  8 |  7 |  6 |  5 |  4 |  3 |  2 |  1 |  0 | -1 |
-x+5 | 10 |  9 |  8 |  7 |  6 |  5 |  4 |  3 |  2 |  1 |  0 |

このとき、y の最小値を達成する直線がどれかを考えます。

   x | -5 | -4 | -3 | -2 | -1 |  0 |  1 |  2 |  3 |  4 |  5 |
2x+3 | ** | ** | ** | ** |    |    |    |    |    |    |    |
 x+1 |    |    |    |(**)| ** | ** | ** |    |    |    |    |
   3 |    |    |    |    |    |    |    |    |    |    |    |
-x+4 |    |    |    |    |    |    |    | ** | ** | ** | ** |
-x+5 |    |    |    |    |    |    |    |    |    |    |    |

特定の直線が最小値となる x は、(高々一つの)区間になっていることがわかります。 また、傾きの降順に見ると、区間は順に並ぶことがわかります。

さらに、傾きが同じ直線が複数ある場合、y-切片が最も小さいものを採用すればよいので、傾きは distinct としてよいです。 これらをまとめて、ざっくり次のようなイメージが湧きます。

  • 「この傾きに対応する最小の y-切片は?」というのを管理する
  • 「この傾きの直線が最小値となる区間はあるか? あるならどこか?」というのを管理する
    • これは、実際には区間の右端だけ持てばよい
    • 区間が空になった傾きは捨ててよい
  • 「この区間で最小値となる直線の傾きはなにか?」というのを管理する
    • 最小値クエリに答えるのに必要になる*1
    • これは一つ上にある「傾き→区間」の逆向きの対応づけ

詳細

左側にある直線 \(y=a_lx+b_l\) と右側にある直線 \(y=a_rx+b_r\) (\(a_l \gt a_r\)) に対して、左側の方が下(同じのも含む)となるのは、 \[a_l x+b_l \le a_r x+b_r \iff x \le \frac{b_r-b_l}{a_l-a_r} \eqqcolon f(a_l, b_l, a_r, b_r)\] となります。\(x\) が整数であれば \(\lfloor\bullet\rfloor\) にしてよいです。参考袋

直線 \(y=a_mx+b_m\) も考えたとき、\(a_l \gt a_m \gt a_r\) に対して \(f(a_l, b_l, a_m, b_m) \ge f(a_m, b_m, a_r, b_r)\) となれば、\(m\) の直線(のみ)が最小となる区間はないことになります。 (雑な記法で書くとして、)\(m \lt r\) となる区間では \(m\ge l\) となっているためです。

というわけで、\( (a, b)\) の追加クエリでは次のようにします。

  • 傾きを distinct にする前処理。
    • \( (a, b')\in S \wedge b' \le b\) であれば何もせず終了。
    • \( (a, b')\in S \wedge b' \gt b\) であれば \( (a, b')\) は消しておく。
  • \(a\) が必要か判定。
    • \(a_l \gt a \gt a_r\) なる \( (a_l, b_l), (a_r, b_r) \in S\) のうち、\(a_l\), \(a_r\) が最も \(a\) に近いものを取る。
      • \(a_l\) か \(a_r\) が取れない場合は \(a\) は必要なので次へ。
      • \(f(a_l, b_l, a, b) \ge f(a, b, a_r, b_r)\) であれば不要なので終了。
        • 上で \( (a, b')\) を消した際はここで終了しないので大丈夫。
  • \(a\) によって不要になった直線の削除。
    • \(a_{l'} \gt a_l \gt a\) を同様に取り、\( (a_l, b_l)\) が不要であれば取り除く。
      • 不要でない \(a_l\) が見つかるまで、取り除くのを繰り返す。
    • \(a \gt a_r \gt a_{r'}\) についても同様に行う。
  • \(a_l\) (if any) と \(a\) が最小となる区間の更新。
    • \(a_l \gt a \gt a_r\) を取り、\(f(a_l, b_l, a, b)\), \(f(a, b, a_r, b_r)\) を計算し、該当の区間を更新する。
      • \(a_l\) が取れない場合、そちらの更新は不要。
      • \(a_r\) が取れない場合、区間の右端は \(\infty\)。

また、\(x\) の最小値クエリでは、\(x\) が含まれる区間における直線 \( (a, b)\) を探し、\(ax+b\) を返せばよいです。

各クエリは map への定数回の処理で対処できるため、\(O(\log(n))\) time です。 より詳細には、クエリ処理前時点における、必要となりうる直線の本数を \(m\) として \(O(\log(m))\) time です。

実装に関して

「傾き→切片」「傾き→区間の右端」「区間の右端→傾き」の 3 つの map を管理することになります。 このとき、「傾きを更新したが区間の右端を更新し忘れた」といったバグを防ぐため、「片方を更新したら他方も合わせて更新する」という処理をする map の wrapper を作っておくとよいかもしれません。

よしなにやって「傾き→(切片, 区間の右端)」「区間の右端→傾き」とすれば空間を若干減らせますが、実装が若干複雑になるかもしれません。

提出コード

  • Library Checker :: Line Add Get Min

感情

CHT で高速化できる DP 面白くて感動しがちだけど、数回見ると飽きがち。

CHT 初めての人は、上に貼った 3 問を考察してみてほしいです。

おわり

にゃんこ。

*1:もう少し考察して、自前の平衡二分木を実装すれば必要なくなるという説もある。

Re: 浮動小数点数の二分探索

rsk0315.hatenablog.com

こういう記事を過去に書いたんですが、どうにもなんだかな〜という気がします。

導入

というのも、64 bit の浮動小数点数を使うとして 90 回もチェックするのはうれしくなさそうです。 (IEEE 754 準拠というのを仮定するとして)double のビット表現は、整数として見ても順序がある程度保たれます。

であれば、整数として見てにぶたんしちゃいたいな〜となります。そうすれば 64 回程度で済むので。

そこで、double としての値とビット表現を並べてみます。

ビット表現
-NaN 1111111111111000000000000000000000000000000000000000000000000000
-inf 1111111111110000000000000000000000000000000000000000000000000000
-1.7976931348623157e308 1111111111101111111111111111111111111111111111111111111111111111
-1.0 1011111111110000000000000000000000000000000000000000000000000000
-2.2250738585072014e-308 1000000000010000000000000000000000000000000000000000000000000000
-5e-324 1000000000000000000000000000000000000000000000000000000000000001
-0.0 1000000000000000000000000000000000000000000000000000000000000000
0.0 0000000000000000000000000000000000000000000000000000000000000000
5e-324 0000000000000000000000000000000000000000000000000000000000000001
2.2250738585072014e-308 0000000000010000000000000000000000000000000000000000000000000000
1.0 0011111111110000000000000000000000000000000000000000000000000000
1.7976931348623157e308 0111111111101111111111111111111111111111111111111111111111111111
inf 0111111111110000000000000000000000000000000000000000000000000000
NaN 0111111111111000000000000000000000000000000000000000000000000000

値が負 (-0.0 含む) であるところに関しては、ビット表現を整数として読んだ順と逆、正 (+0.0 含む) に関しては同じになっていることがわかります。 最上位ビットをよしなにすることも考え、次のような変換を考えます。Rust です。

fn f2u(f: f64) -> u64 {
    let u = f.to_bits();
    if u >> 63 == 1 { !u } else { u | 1 << 63 }
}

値は次のようになります。

f f2u(f) のビット表現
-NaN 0000000000000111111111111111111111111111111111111111111111111111
-inf 0000000000001111111111111111111111111111111111111111111111111111
-1.7976931348623157e308 0000000000010000000000000000000000000000000000000000000000000000
-1.0 0100000000001111111111111111111111111111111111111111111111111111
-2.2250738585072014e-308 0111111111101111111111111111111111111111111111111111111111111111
-5e-324 0111111111111111111111111111111111111111111111111111111111111110
-0.0 0111111111111111111111111111111111111111111111111111111111111111
0.0 1000000000000000000000000000000000000000000000000000000000000000
5e-324 1000000000000000000000000000000000000000000000000000000000000001
2.2250738585072014e-308 1000000000010000000000000000000000000000000000000000000000000000
1.0 1011111111110000000000000000000000000000000000000000000000000000
1.7976931348623157e308 1111111111101111111111111111111111111111111111111111111111111111
inf 1111111111110000000000000000000000000000000000000000000000000000
NaN 1111111111111000000000000000000000000000000000000000000000000000

単調増加になっていることがわかります。それから、NaN が inf より外側にいるのも都合がいいです。

というわけで、これに関してにぶたんすればよいです。

実装例

ABC 026 D を使います。

f2u逆関数 u2f を用意しておきます。

use std::f64::consts::PI;

use proconio::input;

fn main() {
    input! {
        a: f64,
        b: f64,
        c: f64,
    }

    let f = |x: f64| a * x + b * (c * x * PI).sin();

    let res_u = {
        let mut lb = f2u(0.0);
        let mut ub = f2u(1.0e9);
        while ub - lb > 1 {
            let mid_u = lb + (ub - lb) / 2;
            let mid = u2f(mid_u);
            *(if f(mid) < 100.0 { &mut lb } else { &mut ub }) = mid_u;
        }
        ub
    };
    let res = u2f(res_u);

    println!("{}", res);
}

fn f2u(f: f64) -> u64 {
    let u = f.to_bits();
    if u >> 63 == 1 { !u } else { u ^ 1 << 63 }
}

fn u2f(u: u64) -> f64 {
    f64::from_bits(if u >> 63 == 1 { u ^ 1 << 63 } else { !u })
}

atcoder.jp

最上位ビットを持ってきてがちゃがちゃしたいときは、算術シフトを使えば if がいらないんですね。

fn mask(u: u64) -> u64 { ((u as i64 >> 63) as u64 >> 1) | 1 << 63 }

fn f2u(f: f64) -> u64 {
    let u = f.to_bits();
    u ^ mask(u)
}

fn u2f(u: u64) -> f64 { f64::from_bits(u ^ mask(!u)) }

atcoder.jp

人はなぜ

C++ での実装例をよこせと言われそうな気がしますが、C++ は標準で直接 double のビット表現を得る関数がないんですね。bit_cast というのを C++20 まで待つ必要があるらしいです。それから雑にやろうとすると未定義動作にもなります。嫌だねえ。

というので、ちょっと面倒ですが、以下を参考にしつつ書きます。

en.cppreference.com

提出コードのリンクだけ載せておきます。

atcoder.jp

おきもち

定数回ループの方が直感的に思う人も多そうだし、好きな方を使えばいいんじゃないかなとも思います。

それはそれとして、\( (-\infty, +\infty)\) の探索をしたいとき、整数に変換してやれば上から bit を決めていけばよくなって楽かもみたいな気持ちもあります。NaN には注意する必要があるかも。

おわり

ねこにゃんこ。

オフラインの線形時間 LCA

クエリ先読みする場合の LCA を線形時間で解きます。何日か前に TL で話していたら思いついたので。

まずそのときのツイートをいくつか貼ります。

niuez.github.io

クエリ先読みしなくても線形 RMQ を使えば線形時間で解けますが、それはそれです。

定義

念のため。

LCA (lowest common ancestor) problem は次のような問題です。

  • まず根を \(r\) とする \(n\) 頂点の木が与えられて、クエリが \(q\) 個与えられます。
  • \(i\) 個目のクエリ \( (u_i, v_i)\):\(u_i\) と \(v_i\) の両方の祖先である頂点のうち、\(r\) から最も遠い頂点を出力せよ。

クエリ先読みというのは、「クエリを順に読んで一つずつ処理する」のではなく、「クエリを一通り先に読み、まとめてそれらを処理する」という解法の総称です。

オフラインというのはそうした解法を許容する問題設定(あるいはそうしたアルゴリズムも指す)で、そうでないものはオンラインと呼ばれます。競プロにおいてオンラインアルゴリズムが想定とされる(+ オフラインは許容したくない)場合は、インタラクティブ問題として出題されたり、「\(i\) 番目の出力と \(i+1\) 番目の入力の xor を実際の \(i+1\) 番目の入力値とする」などの問題設定で出題されたりします。

解法

まず、適当に DFS するなどしておいて、頂点の番号を pre-order(最初に訪れた順)に振り直しておきます。ここまでは簡単に \(O(n)\) 時間でできます。

      0
     /  \
    1    9
   / \
  2   6
 /|\  |\
3 4 5 7 8

そして、頂点 \(i\) に対応する値 \(a_i\) を管理しながら DFS します。 \(a_i\) は \(1\) で初期化しておき、DFS をして \(i\) から親に帰るとき \(a_i \gets 0\) で更新します。

たとえば、8 に到着した時点の \(a_i\) の値は次のようになります。

      1
     /  \
    1    1
   / \
  0   1
 /|\  |\
0 0 0 0 1

自分より先に訪れた頂点について、自分の祖先であれば 1、そうでなければ 0 になっていることがわかります。

また、表の形式で書くと次のようになります。

i:    0 1 2 3 4 5 6 7 8 9
a[i]: 1 1 0 0 0 0 1 0 1 1

今いる頂点 (8) とそれ以前に訪れた頂点 i について、lca(i, 8) は「(その時点の)a において、j <= i && a[j] == 1 となるような最大の j」、すなわち「i より左(i 含む)に 1 が現れる最も右の位置」となっています。

たとえば、lca(7, 8) == 6lca(4, 8) == 1lca(0, 8) == 0 です(表と木を見比べてみましょう)。lca(9, 8) はこの時点ではまだわかりません。

さて、これらの処理を行うために、以下の処理を \(O(1)\) 時間でできるデータ構造があればよいです。

  • \(a_i \gets 1\) (for all \(i\in[0, n)\)) で初期化。
  • \(a_i \gets 0\) (for given \(i\))で更新。
  • \(\max\,\{j \mid j\le i \wedge a_j = 1\}\) (for given \(i\))を報告。

そういうデータ構造はあります。ABC 228 D 解説 などを参照。

よって、できました。

こうしたデータ構造を持っていなくても、C++ で言うところの std::set などを使うと、\(O(n+q \log(n))\) 時間になります。

実装

judge.yosupo.jp

あんまり定数倍のことは気にしていません。これから考えるかもしれません。

おわり

みんな大好き熨斗袋。もちろんえびちゃんも大好きです。