えびちゃんの日記

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

実数の二分探索・三分探索はループ回数を決め打ちしようね

これ、蟻本とかにも書かれているので常識かなと思ってたんですが、最近よくハマっている人を見かけるので、記事として書いておくことにします。

↓ 追記 ↓

この記事と別の手法についても書きました。競プロの範囲ではどちらの手法でも困らない気もしますが、話として両方知っていても損はしないはず。

rsk0315.hatenablog.com

↑ 追記おわり ↑

まえがき

いつもの

floating-point number を「浮動小数」って呼ぶ or そう呼んで違和感のない人は認識を改めた方がいいと思う。 何が浮動なのかわかってなさそう。

あと「二分探索」は「にぶんたんさく」だと思う。

誤差に関して

誤差があること自体に不満や疑問がある人は、浮動小数点数(や実数表現)について一度勉強するといいかもしれません。 あるいは、必要に応じて任意精度で実数を扱う方法について調べてみると面白いかもしれません。 もちろん、計算量は精度に対して定数ではないので注意しましょう。

本題

探索範囲を \([L, U)\) とし、許容誤差を \(E\) とします。このとき、\(\log_r{\frac{U-L}{E}}\) 回のループが必要になります。 ただし、\(r\) はループ一回ごとの区間の縮小率です。 二分探索なら \(r=2\)、三分探索なら実装によって \(r=3/2\) だったり \(r=\phi=1.618\ldots>3/2\) だったりします。

\(r=\phi\) となる三分探索は、三等分するやつより収束が速かったり、ループ中で計算した値を使い回せたりしてうれしいです。 需要があれば記事を書きますが、とりあえず実装だけ置いておきます。

rsk0315.github.io

実際には、上で計算した値より一回二回くらい多めに見ておくと安心かもしれません。 競技プログラミングの範囲では、\(U-L = 10^{18}\)、\(E=10^{-9}\) とすれば十分でしょうから、\(10^{27}\) についての値を載せておきます。 \[ \begin{aligned} \log_2(10^{27}) &< 90,\\ \log_{\phi}(10^{27}) &< 130,\\ \log_{3/2}(10^{27}) &< 154. \end{aligned} \]

実証

実際に落とせるケースを示さないと納得してくれない人もいるでしょうから、示します。 シンプルな問題でいきましょう。

問題

えびちゃんは \(y\in [0, 10^{18})\) から一つ実数を選びます。 あなたは f(x) と呼び出すことで、\(x \ge y\) の真偽を得ることができます。

あなたには、えびちゃんが選んだ実数を当ててほしいです。 絶対誤差・相対誤差の大きくない方が \(10^{-9}\) 以下であれば許容されます。

解法?

あなたは次のようなコードを書きました(書くでしょう?)。

int main() {
  long double lo = 0.0L, hi = 1e18L;
  while (hi-lo > 1e-9L) {
    long double mid = 0.5L * (lo+hi);
    (f(mid)? hi: lo) = mid;
  }
  printf("%.21Lf\n", lo);
}

落とすよ

落とします。 えびちゃんは次のように f を定義します。

bool f(long double x) {
  return x >= 1e17L;
}

Successful hacking attempt (tl)

理由の説明

精度の関係で、ループの途中からは次のような状態になっています。

lo == 99999999999999999.9921875L
hi == 100000000000000000.0L

これらの差は、0.0078125L ですから、< 1e-9L を満たしません。

long double において、この lo より大きい最小の値は hi なので、ループ継続条件は偽になりません。 大抵の場合は、こうなると無限ループになります*1

lo より大きい最小の値を知りたければ、次のようにするとよいです。

#include <cmath>
#include <cstdio>
#include <limits>

int main() {
  long double lo = 99999999999999999.9921875L;
  printf("%.21Lf\n", std::nextafter(lo, std::numeric_limits<long double>::infinity()));
}

こんな感じです。

ループ継続条件に相対誤差も考慮して、

while ((hi-lo) / std::max(1.0L, lo) > 1e-9L)

とかするといいかもなんですが、どうでしょう。割り算あんまりしたくないしね。

今回は、煩雑になるのを避けるために、探索したい値をえびちゃんが選んだ設定にしましたが、「何らかの判定関数に関して、境界が \(y\) であった」というのと同じことです。

三分探索でも似たように落とせると思います。

思ったんだけどさ

下に凸な \(f(x)\) の最小値を求めてねって問題を考えましょう。

\(x^\ast = \arg\,\min_x f(x)\) とします。 で、long double において \(x^\ast\) より小さい最大の数を \(x_-\)、それより大きい最小の数を \(x_+\) とします。 そして、\(f\) が非常に厄介で、\(f(x_+)-f(x^\ast)\) や \(f(x_-)-f(x^\ast)\) がとっても大きいとします。

...というの、こわれると思うんですが、作った問題がうっかりこうなってないか気をつけなきゃですね。

それにちょっと関連して、\(x\) の範囲を十分小さめに絞れたとしても、その区間での \(y\) の値域が十分小さいとは限らないので注意が必要そうです。 このことは、範囲の大きさでループするにしても、ループ回数を決め打ちするにしても、気をつける必要がありそうです。

*1:言い切っていないのは、こういうループの動作は未定義だからです。未定義動作についての記事