えびちゃんの日記

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

ならし計算量のよくある誤解について

ならし計算量 (amortized complexity) というのがありますが、競プロ初心者層には誤解されがちな概念なので、誤解を解くためのことを書きたくなりました。

以下で「よくある誤解」と称して架空の人物がしゃべりますが、複数回見てきた事例であり、特定個人に対する揶揄等ではありません。えびちゃん自身の事例である場合もあります*1

多少知っている人へ:ポテンシャル解析などの話を直接はしないので、そういうのは別の記事を見に行ってほしいです。文献の紹介くらいはあります。

導入

データ構造の計算量などで、次のような記述を見かけることがよくあります。

  • vector へのアクセス a[i] は、最悪 \(O(1)\) 時間
  • set での二分探索 a.lower_bound(x) は、最悪 \(O(\log(n))\) 時間

これらは、「一番悪いときでも \(O(1)\) 時間や \(O(\log(n))\) 時間だから安心だな〜」という気持ちで読める人が多そうです。

一方で、次のような記述もよくあります。

  • vector の末尾への追加 a.push_back(x) は、ならし \(O(1)\) 時間
    • 償却 \(O(1)\) 時間という言い方もあり、これらは同じ意味
    • 同様に、ならし計算量を償却計算量とも言う
  • union_find での操作 a.unite(i, j)a.find(i) は、ならし \(O(\alpha(n))\) 時間

これらは、誤解して「ならしってなに? やばいときは \(O(1)\) 時間や \(O(\alpha(n))\) 時間にならなくて、知らないうちに TLE になったりしそう?で不安だな〜」という気持ちになる初心者が多い印象があります。

実際には、少なくとも競プロ文脈で言えば、「ならしで \(O(1)\) 時間や \(O(\alpha(n))\) 時間だから安心だな〜」となれるくらいのうれしい保証なので、その説明をします。

説明

お気持ち

データ構造に対して \(n\) 回操作し、\(n\) 回の合計の計算量が \(f(n)\) であったとき、一回あたりの計算量は \(f(n)/n\) と見なせます。 そこで、実際の操作一回一回の計算量は考えず、この \(f(n)/n\) について考えてみることにします。

たとえば \(f(n) = O(n)\) であれば、一回あたり \(O(1)\) 時間ということになります。

これをもう少しちゃんと保証したものが、ならし計算量と呼ばれるものです。

定義(ゆるふわ)

あるデータ構造に対して行える操作(vector に対する push_back とか、stack に対する pushpop とか)を考えます。それらをどういう順番で \(n\) 回処理した場合でも、合計が \(f(n)\) 時間以下になるとします。このとき、それらの操作は「ならし \(O(f(n)/n)\) 時間」であると言います。

たとえば、「union_find に対する操作はならし \(O(\alpha(n))\) 時間」というのは、「union_find に対して \(n\) 回操作をすると \(O(n\cdot \alpha(n))\) 時間」のような意味合いです。

定義(かっちりめ)

データ構造に対する任意の操作列 \( (op_1, op_2, \dots, op_n)\) を考えます(\(op_1\) が .push_back(x) とか、そういう感じです)。 このとき、\(op_i\) (\(1\le i\le n\)) に対して、実際の計算量を \(t(op_i)\) とします。 さらに、\(a(op_i)\) が存在して、各 \(1\le k\le n\) に対して以下が成り立つとします。 \[ \sum_{i=1}^k t(op_i) \le \sum_{i=1}^k a(op_i). \] このとき、\(a(op_i)\) を(\(op_i\) の)ならし計算量と呼びます。

\(t(op_i) \le a(op_i)\) (\(1\le i\le n\)) が成り立つ必要はなく、\(i\) 回目以前の操作の合計さえ抑えられればよいのがミソです。

\(n\) 回の操作全体の計算量が \(\sum_{i=1}^n a(op_i)\) で抑えられるため、たとえば \(op_i\) が「ならし \(O(\log(n))\) 時間」と保証されている操作であれば、全体で \(O(n\log(n))\) 時間とわかったりするわけです。

(補足) この \(a(op_i)\) をうまいこと求めるのが、各問題でならし計算量を求めるために必要なことであり、ならし解析 (amortized analysis) などと呼ばれているものです。解析の具体例だけを見せられると天才に見えますが、典型的な手法がいくつかあります。

やや細かい話

実際には、\(n\) が「操作列の長さ」なのか「その時点のデータ構造のサイズ」なのか「それ以前のデータ構造のサイズの最大値」なのかがぼかして話されたりしがちな気もします。

「空の状態から始めて、操作ごとに増える要素の個数が高々 1 つ」であるような場合には、データ構造のサイズが操作列の長さで上から抑えられるので、問題なさそうな気もします。

そうでない場合は少し気にする必要があるかもしれません。 「サイズ \(n\) で初期化して \(m\) 回操作する」ような場合には、初期化の計算量 \(t(op_0)\) も合わせて考えるかもしれません。 あんまり初期化のことを書いている文献はないかもしれません。

\(\Theta\) とか \(\Omega\) とかの話

ならし計算量の定義だと上からしか抑えてないけど、ならしで \(\Theta\) や \(\Omega\) ってうまく定義されてる?

「ある操作列が存在して」か「任意の操作列に対して」かは文脈にもよりそうだけど、下から抑えれば定義できる? できそう。

\(O\) 以外を知らない人向け → やや長めの記事

よくある例

可変長配列

vector.push_back のようなものです。ざっくり以下のようなアルゴリズムです。

  • 初め、長さ \(1\) の領域を確保する
  • 確保した領域に新しく要素を置けるなら置く(\(O(1)\) 時間)
  • 確保した領域が足りなくなったら、長さを倍にしてから置く(領域の長さを \(k\) として \(\Theta(k)\) 時間)

長さ \(n\) にするための計算量は?となり、「\(\Theta(n)\) 時間の操作があって \(n\) 回だから \(\Theta(n^2)\) 時間か...?」とか「\(\Theta(n)\) 時間の操作は \(\floor{\log_2(n)}\) 回起きるから \(\Theta(n\log(n))\) 時間か...? \(\log\) は定数()だから一回あたり \(O(1)\) 時間と言われているだけなのか...?」となったりしますが、そうではないです。

長さ \(k\) は毎回固定ではなく \(1 + 2 + 4 + \dots + 2^{\floor{\log_2(n)}}\) のようになっているため、(等比数列の和であることに注意するなどして)全体で \(\Theta(n)\) 時間であることがわかります。

キュー

二つのスタックでキューを作ることができ、two-stack queue などと呼ばれます*2

下の方にあるスライドにも載っているので、詳細は割愛します。似た方法で deque も作れます*3

カウンタ

\(0\) から始めて \(n\) まで \(1\) ずつ増やすことを考えます。このとき、各インクリメントに対して、数字が変化する桁はならし \(O(1)\) 個であることが示せます。 なので、変わった桁のみをうまく管理することで、たとえば以下のようなクエリをならし \(O(1)\) 時間で処理できます。

  • 値を \(1\) 増やす
  • 桁和を出力する
    • あるいは桁の二乗和とかを出力する

ゆるふわな話

ならしのイメージの話です。

たとえば、「1 日に使った金額が 100 円を超えちゃいけませんよ」と「1 週間に使った金額が 700 円を超えちゃいけませんよ」という制限を考えます。 後者は「平日は 1 円も使わずに土日に 350 円ずつ使う」といった行動が許されるのに対し、前者ではそうしたことが許されません。

この前者が最悪計算量(各クエリの実際の計算量)に相当し、後者がならし計算量(全体で見たときに一回あたりと見なせる計算量)に相当するつもりです。

より実態に即しているイメージをフォロヮにもらったので書きます。前者は「毎日 100 円ずつあげるけど、使わなかったぶんは没収ね」で、後者は「毎日 100 円ずつあげるから、使っても貯金してもいいよ(借金はだめ、利子はなし)」です。計算量が小さい操作をしたのが貯金を作ったことに相当します。

ならしがうれしい話

ならし計算量が解析の面でうれしい状況としては、

  • 一回一回の解析は難しいが、全体で見れば解析が簡単
  • 重いクエリもあるが、全体で見れば十分軽い

のようなものがあります。それはそれとして、「全体としての計算量さえ保証すればいいときは、ならしだけ保証するように設計する」という方針がうれしいケースもあります。

たとえば、\(n\) 個の操作の計算量が \( (1, 1, \dots, 1, n)\) で合計 \(2n-1\) になっているデータ構造があったとします。 これはならし \(O(1)\) 時間ですが、最悪 \(O(1)\) 時間ではありません(最後の一回が \(n\) になっているので)。

これをなんとか工夫することで \( (10, 10, \dots, 10, 10)\) にできたとします。これは最悪 \(O(1)\) 時間を保証できていますが、合計としては \(10n\) に増えてしまっています。

このとき、定数倍が悪くなっている上、全体のオーダーが改善されているわけでもなく、実装の工夫も必要になってしまい、大変づくしです。

もちろん、最悪計算量を保証するために必ずしも定数倍がハチャメチャに悪化するわけではないので状況によりけりですが、ゆるい制限でいいならそれに応じてうまくやりましょうということです。

ならしだと困る話

競プロだとあんまりないと思います。

たとえば、Web サイトなどで「\(i\) 番目のアクセスに対しては、データ構造に対して処理 \(i\) を行う」という処理があったとします。 このとき、ならし計算量しか保証されていないと、一部のユーザが「めちゃくちゃ時間かかった」というような体験をすることになり、あまりうれしくなさそうです。

競プロでも困る例をフォロヮにもらったので書きます。 「ならし○○時間」の操作を普通にするだけなら問題ないと思いますが、rollback が可能なデータ構造に書き換えようとするときに、元が「ならし」のデータ構造だとつらいことがあります。

先のカウンタの例で言えば、999...9 から(大きい計算量で)1000...0 に更新したあと、rollback をして(大きい計算量で)999...9 に戻り、また(大きい計算量で)1000...0 に更新して、... というのを繰り返すと、毎回大きい計算量がかかるようになりますね。

とはいえ、「ならし○○時間」と保証されているデータ構造を、保証されている通りに使えば問題ないことには変わりない気はします。

ならし計算量あるあるデータ構造

C++/STL にもいろいろありますね。

  • vector.push_back(x) はならし \(O(1)\) 時間
    • このせいで、vector を内部で使うデータ構造は大抵ならしになりがち
  • set.erase(it) はならし \(O(1)\) 時間
    • コストを .insert(x) に押しつけることで、.erase(it) の方はならせそう?
    • それはそれとして、赤黒木を平衡させる操作はならし \(O(1)\) 時間
      • どこに挿入・削除するかは、ポインタなどを持っていないと \(O(1)\) で決められないことに注意
      • itイテレータ

STL ではないものの競プロでよくあるデータ構造にもいろいろあります。

  • union_find.unite(i, j).find(i) は(適切に実装すれば)ならし \(O(\alpha(n))\) 時間
    • 実際には、サイズ \(n\) のそれに \(m\) 回操作すると \(\Theta(n+m\,\alpha(m, n))\) 時間で、\(O(\alpha(n))\) 時間よりも厳しく抑えられる
    • \(\alpha(m, n)\) については → そういう記事
  • foldable_queue.fold().push(x) はならし \(O(1)\) 時間
    • いわゆる SWAG(を勝手にそう呼んでいます)
      • queue への push/pop と、全体のモノイド積 (fold) を処理するデータ構造
      • push と fold は最悪 \(O(1)\) 時間
    • two-stack queue を応用して実装する(ので deque でも可能)
  • decremental_neighbor_query.less_than(i).greater_equal(i) などは、ならし \(O(1)\) 時間
    • \(\{1, 2, \dots, n\}\) で初期化した集合に対し、削除操作と、隣の要素の検索ができる
  • skew_heap.meld(q) はならし \(O(\log(n))\) 時間
    • ヒープ(優先度つきキュー)二つに対して、それらをくっつけて一つのヒープにする操作を meld と呼ぶ
  • fibonacci_heap.prioritize(it, k) はならし \(O(1)\) 時間
    • ヒープに対して、ある要素の優先度を高める操作(取り出されやすくする)を prioritize と勝手に呼んでいます
      • 文献によっては decrease_key とか言うかも
    • この操作を \(O(1)\) 時間でできることにより、Dijkstra 法の計算量を \(O(|E|\log(|V|))\) から \(O(|E|+|V|\log(|V|))\) に改善できる

まとめ

ならし計算量は、操作列を全体として見たときの計算量を保証してくれるため、競プロ文脈では「最悪○○時間」と同程度にえらい保証と見なしてよいです。 「入力によっては○○時間になってくれず TLE(平均計算量)」とか「乱数の値によっては○○時間になってくれず TLE(期待計算量)」とかとは事情が異なります。

参考文献

  • cs166.1266
    • 図とかがたのしい感じの講義スライド
    • 他にもいろいろなトピックがある
  • 熨斗袋の記事
    • ならし以外にも、期待計算量や平均計算量についても載っている

おわり

にゃんこ。

*1:そういう意味では特定個人ですが、自分なので許されます。

*2:文法としては、two-week vacation とかのように、ハイフンで数詞とつないで名詞が単数形になる形容詞のやつです。two-stacks queue や two-stack-queue ではないです。

*3:概要:片方が空になったときに \(k\) 個移すのではなく \(\ceil{k/2}\) 個移します。

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

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

動機について

やや前に 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 できそうな例を出して初心者を脅かしたかったんですが、いいのが思いつきませんでしたね。

sqrt

\(0\le n\le 10^{18}\) で \(\floor{\sqrt{n}}\) を求めてね、と言われて、std::sqrt(n) と書くとだめです。 sqrt が悪いというよりは、double を受け取る関数に long を渡した際に暗黙に型変換が起きているのが原因なのですが。

たとえば \(n = 2^{54}-1\) のとき、double では \(2^{54}\) に丸められ、\(2^{27}\) が返されたりします。実際には \(\floor{\sqrt{n}} = 2^{27}-1\) なので嘘になりますね。 平方数より少し小さい整数で、上に丸められてその平方数以上になるケースであれば、なんでも hack できると思います。

指数探索なり二分探索なり Newton 法なり整数でやる方針にするか、double で返ってきた値の周辺を線形探索する方針がよさそうです。

おまけ

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

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

敢えて作ったような例ですが、「\(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 の略です。どこ発祥の呼び方なんだろう。

上記のクエリを高速に処理するデータ構造を指して CHT と言っている人や、特定の DP を上記のクエリに帰着するテクを指して CHT と言っている人がいる気がします。えびちゃん的には後者の気持ちなのですが、ここでは便宜上前者の意味で言っていそうです。

アイディア

たとえば、\(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

その他

github.com

C++ での実装です。multiset を使ったり、mutable を利用していたりイテレータをがちゃがちゃやっていて少し怖いですが、実装量はとても短いです。

kactl は、KTH (Kungliga Tekniska högskolan; Royal Institute of Technology) のチームの ICPC 用ライブラリで、内容が充実しています。

感情

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 には注意する必要があるかも。

おわり

ねこにゃんこ。