浮動小数点型の二分探索についてはいろいろと議論があります。
背景
下記のような形式の問題を考えます。
実数に関する述語 $P\colon \R\to\{0, 1\}$ と $x_L, x_R\in\R$ が与えられる。
ある $x^{\ast}\in[x_L\lldot x_R]$ が存在し、任意の $x\in[x_L\lldot x_R]$ に対して $P(x) = 1 \iff x\le x^{\ast}$ が成り立つことが知られているため、その $x^{\ast}$ を求めよ。
ただし、$x^{\ast}$ との絶対誤差または相対誤差が $10^{-9}$ 以下である値を出力すれば正答と見なす。 また、$x_R-x_L \le 10^{18}$ とする。
無邪気な競プロ er*1は、浮動小数点型を用いて、次のような二分探索を実装しがちです。 言語による差異は特にないですが、今回は Rust で書いてみます。
fn float_bisect_1(mut xl: f64, mut xr: f64, p: impl Fn(f64) -> bool) -> f64 { while xr - xl >= 1.0e-9 { let xm = 0.5 * (xl + xr); *(if p(xm) { &mut xl } else { &mut xr }) = xm; } xl }
たとえば $(x_L, x_R) = (0, 10)$ かつ $P(x) = 1 \iff x^2\le 3$ であれば
float_bisect_1(0.0, 10.0, |x| x * x <= 3.0)
で 1.732050807331688702106475830078125 となり、$\sqrt3$ に近くてよさそうな見た目をしています。しかし、次のような例ではうまくいきません。
float_bisect_1(0.0, 2.0e7, |x| x * x <= 1.0e14)
二分探索の過程で (xl, xr) = (10000000.0, 10000000.00000000186264514923095703125) となったとき、丸め誤差の関係で 0.5 * (xl + xr) == xl となってしまいます。
すなわち xm == xl となるため、結果として無限ループになってしまいます*2。
蟻本の時代からある定跡としては、定数回のループを行うことで無限ループを避けるというものです。$x_R-x_L = 10^{18}$ と許容誤差の $10^{-9}$ を考慮して、$\log_2(10^{18}\cdot 10^9) \lt 90$ 回程度の更新で十分そうです。
fn float_bisect_2(mut xl: f64, mut xr: f64, p: impl Fn(f64) -> bool) -> f64 { for _ in 0..90 { let xm = 0.5 * (xl + xr); *(if p(xm) { &mut xl } else { &mut xr }) = xm; } xl }
こちらは次のようにして期待通りの結果が得られます。
float_bisect_2(0.0, 10.0, |x| x * x <= 3.0) // 1.732050807568877193176604123436845839023590087890625 float_bisect_2(0.0, 2.0e7, |x| x * x <= 1.0e14) // 10000000.0
ではこれでめでたしめでたしでしょうか？というのが今回つっつく重箱の隅です。
ここまでの方法においては、停止性についてしか議論しておらず、部分正当性についての議論をしていません。
すなわち、ループが終了したときに得られる値が $x^{\ast}$ と十分に近いことが言えていません。
特に、浮動小数点型で実装した p が $P(x) = 1 \iff x\le x^{\ast}$ の性質を持っていることを暗に仮定してしまいがちです。
本題
安易にやるとだめになってしまうケースを考えていきます。
単純な撃墜ケース
まずは単純なところからです。 $(x_L, x_R) = (1, 1600)$ かつ $$ P(x) = 1 \iff \bigl(x\le 1000\bigr) \wedge \left(\tfrac1x\cdot x\ge 1\right) $$ とします。あからさまな悪意がむんむんです。これは単に $x\le 1000$ と同値ですが、特に式変形などをせずにそのまま実装すると次のようになります。
let p = |x| x <= 1000.0 && 1.0 / x * x >= 1.0; float_bisect_2(1.0, 1600.0, p); // 775.5155233629557187668979167938232421875
( ﾟдﾟ)ﾎﾟｶｰﾝ
// 775.5155233629557187668979167938232421875
(つд⊂)ｺﾞｼｺﾞｼ
775.5155233629557187668979167938232421875
これはだめですね。
f64 においては
1.0 / 800.5 * 800.5 == 0.99999999999999988897769753748434595763683319091796875
が成り立つのが重要ポイントです。二分探索の過程で調べることになった他のいくつかの値についても、同様に 1.0 / x * x < 1.0 が成り立ちました。
浮動小数点型の撃墜ケースを作る際には、$x/x$ や $x-x$ のような「本来であれば打ち消し合う項」が狙い目のひとつです。
それっぽい例
単純な例（悪意があからさますぎて実践では出てこなさそうな例）はできたので、あとはカモフラージュするだけです。
1.0 / 49.0 * 49.0 も上記と同じ値になる*3ことを踏まえて、次のような問題を考えてみましょう。
多少無理があるかもしれませんが許してもらいます。
多項式 $p(x)$, $q(x)$, $r(x)$, $s(x)$ が与えられる。 $\frac{p(x)}{q(x)}\cdot r(x)\ge s(x)$ なる $x$ の最大値 $x^{\ast}$ を求めよ。
ただし、$x^{\ast}$ との絶対誤差または相対誤差が $10^{-9}$ 以下であれば正答と見なす。
制約
- $p$, $q$, $r$, $s$ は高々 $3$ 次式で、係数の絶対値は $300$ 以下の整数である。
- $0\le x^{\ast}\le 4$ が成り立つ。
- $x\le x^{\ast}$ ならば $\frac{p(x)}{q(x)}\cdot r(x)\ge s(x)$ が成り立つ。
そして、次の入力を考えます。 $$ \begin{aligned} p(x) &= 2x^2-5x+4, \\ q(x) &= 49x^2-98x+98, \\ r(x) &= 147x^2-294x+196, \\ s(x) &= 4x^3-12x^2+11x-2. \end{aligned} $$
さて、これらの多項式について次が成り立ちます。
| $x$ | $p(x)$ | $q(x)$ | $r(x)$ | $s(x)$ |
|---|---|---|---|---|
| $1$ | $1$ | $49$ | $49$ | $1$ |
| $2$ | $2$ | $98$ | $196$ | $4$ |
そのため、f64 で計算すると次のようになります。
let p = |x| (2.0 * x - 5.0) * x + 4.0; let q = |x| (49.0 * x - 98.0) * x + 98.0; let r = |x| (147.0 * x - 294.0) * x + 196.0; let s = |x| ((4.0 * x - 12.0) * x + 11.0) * x - 2.0; (p(1.0) / q(1.0) * r(1.0), s(1.0)) // (0.99999999999999988897769753748434595763683319091796875, 1.0) (p(2.0) / q(2.0) * r(2.0), s(2.0)) // (3.999999999999999555910790149937383830547332763671875, 4.0)
そのため、$[0\lldot 4]$ の範囲で二分探索を開始すると、$x^{\ast}\le 1$ だと思い込んでしまうことになります。 さて、これは正しいでしょうか？
あらあら、正しくなさそうです。 $$ \begin{aligned} &\phantom{{}={}} \frac{p(x)}{q(x)}\cdot r(x)-s(x) \\ &= \frac{2x^2-5x+4}{49x^2-98x+98} \cdot (147x^2-294x+196) - (4x^3-12x^2+11x-2) \\ &= \frac{-2(x-2)(x-1)^2\, (2x^2-5x+5)}{x^2-2x+2} \end{aligned} $$ であり、$2x^2-5x+5=0$ は実数解を持たないことに注意すると、赤実線と黒破線の共有点は $x=1$ と $x=2$ のみであることがわかります。分母が常に正であることにも気をつけておくといいかも。
ということで、実際には $x^{\ast} = 2$ です。1.0 と 2.0 のとき以外にも計算誤差の影響により、大小比較の結果が実数での比較とは異なることが数回あり、最終的には次の値が返ってきました。
let p = |x| (2.0 * x - 5.0) * x + 4.0; let q = |x| (49.0 * x - 98.0) * x + 98.0; let r = |x| (147.0 * x - 294.0) * x + 196.0; let s = |x| ((4.0 * x - 12.0) * x + 11.0) * x - 2.0; let x = float_bisect_2(0.0, 4.0, |x| p(x) / q(x) * r(x) >= s(x)); float_bisect_2(0.0, 4.0, |x| p(x) / q(x) * r(x) >= s(x)) // 0.999999999999999555910790149937383830547332763671875
多項式の計算のところで Horner’s scheme を使わずに次のようにした場合も、同種の誤差でさらに離れた値が返ってきました。
let p = |x| 2.0 * x * x - 5.0 * x + 4.0; let q = |x| 49.0 * x * x - 98.0 * x + 98.0; let r = |x| 147.0 * x * x - 294.0 * x + 196.0; let s = |x| 4.0 * x * x * x - 12.0 * x * x + 11.0 * x - 2.0; float_bisect_2(0.0, 4.0, |x| p(x) / q(x) * r(x) >= s(x)) // 0.9999999999999993338661852249060757458209991455078125
note: 次の入力の場合、$1-6.02\cdot 10^{-9}$ 程度の値が返ってきました。すなわち、$1$ との誤差ですら $10^{-9}$ より大きいです。$q(x) = 49p(x)$ でつまらないかも？ 探せばもっといろいろありそうです。 $$ \begin{aligned} p(x) &= 2x^2-5x+4, \\ q(x) &= 98x^2-245x+196, \\ r(x) &= 98x^2-147x+98, \\ s(x) &= 2x^3-6x^2+7x-2. \end{aligned} $$
他の例
区分線形関数などもよい例かもしれません。たとえば $$ f(x) = \begin{cases} \tfrac{3}{187}(x-13) - 2, & \text{if}\; x\le 200, \\ \tfrac{3}{187}(200-13) - 2, & \text{if}\; 200 \le x\le 250, \\ \tfrac{2}{187}(x-250) + 1, & \text{if}\; x\ge 250 \\ \end{cases} $$ として、$f(x) \le 1$ なる $x$ の最大値を求めてみましょう。範囲は $[0\lldot 400]$ としてみます。
let f = |x| { if x <= 200.0 { 3.0 / 187.0 * (x - 13.0) - 2.0 } else if x <= 250.0 { 3.0 / 187.0 * (200.0 - 13.0) - 2.0 } else { 2.0 / 187.0 * (x - 250.0) + 1.0 } }; float_bisect_2(0.0, 400.0, |x| f(x) <= 1.0); // 199.999999999999971578290569595992565155029296875
3.0 / 187.0 * 187.0 == 3.000000000000000444089209850062616169452667236328125
が成り立つことなどから変な感じになっていき、$200$ より少し小さい値を出力してしまいます。
なにかの部分問題として「こういう区分線形関数を考えて、それがこの値以下である範囲を求めればよいですね〜」となったときに、油断してやらかしてしまいそうな気配を感じます。
まとめ
二分探索を行う際の決定問題パートを素朴にやると失敗してしまう例を挙げました。 二分探索が悪いというよりは、決定問題への帰着を素朴に未証明でやるのが悪いという話ではあります。
単純な例では（いわゆる）eps を足し引きする手法が使えるかもしれませんが、それを行う場合でも正当性の証明はきちんとやる必要があるでしょう。
コンテスタント目線では AC するまで適当にがちゃがちゃやる手法が一応可能ではありますが、作問者側が気づかずに踏んでしまうと大変なことになってしまいそうです。 よくある「想定解との誤差が何々以下のとき AC とします」でわくわくしそうなやつですね。
おまけ
f64 は 64-bit なので、述語の判定を 64 回程度で済ませられるのでは？というトピックもあります。
判定の際にどの値を使うかを決めるパートが異なるだけで、今回の問題と同様の部分に気をつける必要があります。
あとがき
こうした部分に触れられることってあまりないなというのは昔から思っていましたが、1.0 / x * x >= 1.0 のような単純なものしか考えたことがなかったので、今回改めて考えてみました。
Rump の例題もそうですが（そうだと思っていますが）、本質パートを考えたらあとはそれっぽくカモフラージュしていくパートになり、そこが楽しかったり大変だったりしますね。
気づいたらまた浮動小数点型の記事を書いており、えびちゃんはもしかして浮動小数点型が好きなのですか？
おわり
おわりです。
