Sqrt Inequality の誤差評価

atcoder.jp

å°Žå…¥

これの long double 解法についての話です。下記のような解法です。

#include <cmath>
#include <cstdio>

constexpr long double eps = /* ??? */;

bool ok(long double a, long double b, long double c) {
  return std::sqrt(a) + std::sqrt(b) + eps < std::sqrt(c);
}

int main() {
  long double a, b, c;
  scanf("%Lf %Lf %Lf", &a, &b, &c);
  puts(ok(a, b, c) ? "Yes" : "No");
}

ただし、long double の仮数部は 64-bit であるとします（AtCoder 環境に準拠）。 公式解説 では 1.0E-14 に設定し、

ここで $\varepsilon = 10^{-14}$ の値を大きくしすぎたり小さくしすぎたりすると WA になってしまいます。 なぜ $\varepsilon = 10^{-14}$ であるとうまくいくのかは説明できるのですが、より簡単で安全な方法は、全て整数でやってしまうことです。

と書かれています（1.0E-14 は double のリテラルであり、long double の 1.0E-14L より若干精度が劣りますが、ここでは問題になりません）。実際、

eps = 8.8817841970012523238e-16L;

では WA になり、

eps = 8.8817841970012523239e-16L;

では AC になります。また、

eps = 1.5099033134902128948e-14L;

では AC になり、

eps = 1.5099033134902128949e-14L;

では WA になります。

ざっくり 4 桁で見ると、$8.882\times 10^{-16} \le \text{\texttt{eps}} \le 1.509\times 10^{-14}$ の範囲であれば AC になります。

数値のざっくりとした出処は、$2^{\floor{\log_2(\sqrt{10^9})}-64} = 2^{-50}$ と $$ \underbrace{(\ceil{\tfrac12(10^9)^{-3/2}\cdot 2^{(64-1)-\floor{\log_2(\sqrt{10^9})}}}-\tfrac12)}_{8.5} \times \underbrace{2^{\floor{\log_2(\sqrt{10^9})}-(64-1)}\vphantom{\tfrac12}}_{2^{-49}} $$ です（これらに真に含まれる範囲）。

これをちゃんと考えていきましょうというのが今回のお話です。なお、後述しますが、上記の eps = 1.5099033134902128948e-14L による解法は hack することが可能です。

定義・前提知識

$\gdef\sqrtf{\text{\texttt{sqrt}}}$ $\gdef\roundf{\operatorname{round}}$

浮動小数点数型についての基本的な知識はあるものとします（表現の仕方、丸めによって誤差が生じることなど）。ない人のための記事はそのうち書く予定ではあります。

実数 $x$ に対し、（対象としている浮動小数点数型で）表せる最も近い数値*1を $\roundf(x)$ および $\hat x$ で表します（$x$ の部分が長いときには前者の記法を好んで使い、主に一文字の変数のときには後者の記法を好んで使うことにします）。 浮動小数点数型の数値 $x$, $y$ に対して、$x\oplus y = \roundf(x+y)$ と $\sqrtf(x) = \roundf(\sqrt{x})$ と書くことにします。

和および平方根の計算では、正確に計算した値を正確に丸めた値を返すことが定められているので、誤差は 0.5 ULP です。すなわち、正確な値が $2^e$ 以上 $2^{e+1}$ 未満の範囲にあるとき、誤差の絶対値は $2^{e-p}$ 以下です。 ここで、AtCoder 環境における long double においては $p = 64$ です。double や一部環境の long double では $p = 53$、別の環境の long double では $p=113$ だったりします。

また、正の実数 $x$, $y$, $x+y$ が正規化数の範囲に入るとき、$\roundf(x+y)$ と $(\hat x\oplus\hat y)$ の差の絶対値は 1 ULP 以下であることが示せます。すなわち、「丸めた値同士を足してから丸めた値」と「無限精度で足したものを一度だけ丸めた値」の差が 1 ULP で抑えられるということです*2。直感的には、$\roundf(x+y)=(\hat x\oplus\hat y)$ であるか、浮動小数点数型で表せるうちで隣り合う 2 数になります。

$\gdef\ulp{\operatorname{ulp}}$

実数 $x$ に対し、$\ulp(x) = 2^{\floor{\log_2(|x|)}+1-p}$ とします。$2^i\le x\lt2^{i+1}$ のとき、$\ulp(x) = 2^{i+1-p}$ となります。 また、 $$ \ulp(x)\cdot 2^{p-1}\le x\lt \ulp(x)\cdot2^p $$ や $$ \ulp(x)\cdot 2^{p-1}\le \hat x\le \ulp(x)\cdot2^p $$ が成り立ちます。

方針

誤差のない世界で計算した結果と、誤差のある世界で計算した結果とで、大小関係が常に同じだと楽です（誤差を無視して計算して答えればよいので）。 しかし、実際には次のようなケースが存在します。

$$ \sqrt{a} + \sqrt{b} = \sqrt{c} \wedge \sqrtf(a) \oplus \sqrtf(b) \lt \sqrtf(c). $$

そこで、適切な $\varepsilon$ を足すことで、上記の「大小関係が常に同じ」を成り立たせる方法が知られています（蟻本 p. 228 など）。 $$ \sqrt{a} + \sqrt{b} \lt \sqrt{c} \iff (\sqrtf(a) \oplus \sqrtf(b)) \oplus \varepsilon \lt \sqrtf(c). $$ ここで、平方根を求めた際やそれらの和を求めた際だけでなく、$\varepsilon$ を足した際にも丸め誤差が生じることに注意しておきましょう。

非常に嘆かわしいですが、この $\varepsilon$ をまともに考えて設定している人は、中級者以上でも非常に少ない印象があります*3。次のようなことを考える必要があるでしょう。

まず、$\sqrt{a} + \sqrt{b} \ge \sqrt{c}$ のケースにおいて、$\sqrtf(a) \oplus \sqrtf(b)$ が $\sqrtf(c)$ をどの程度下回るかを考える必要があります。 「$\sqrt{a}+\sqrt{b}\ge\sqrt{c}$ のケースで No と答えるために、どの程度 $\varepsilon$ を大きくする必要があるか？」ということです。

また、$\sqrt{a} + \sqrt{b} \lt \sqrt{c}$ のケースにおいて、$\sqrtf(a) \oplus \sqrtf(b)$ が $\sqrtf(c)$ にどの程度近づくかを考える必要があります。 「どの程度なら $\varepsilon$ を大きくしても、$\sqrt{a}+\sqrt{b}\lt\sqrt{c}$ のケースで Yes と答えるままでいられるか？」ということです。

未知の $\varepsilon_L$, $\varepsilon_U$ に対して $\varepsilon_L \le \varepsilon \le \varepsilon_U$ の範囲で AC できるとしたときに、$\varepsilon_L$ の上界 $\varepsilon_L'$ を求めるのが前者で、$\varepsilon_U$ の下界 $\varepsilon_U'$ を求めるのが後者です。$\varepsilon_L'\le\varepsilon_U'$ となるように求められれば、その範囲の $\varepsilon$ を使うことで十分となります。

逆に、そうした区間 $[\varepsilon_L\ldots\varepsilon_U]$ が空の場合はこの解法では正答できないことになります。

考察

一般性を失わず、$a\le b$ とします。また、制約の範囲において $a$, $b$, $c$ は注目している浮動小数点数型で正確に表せるので、そのことも仮定します。

$\varepsilon_L'$ の上界

$\sqrt{a}+\sqrt{b}\ge\sqrt{c}$ とし、$(\sqrtf(a)\oplus\sqrtf(b))-\sqrtf(c)$ を下から抑えます。 丸めを一度しか行わない場合の値と隣り合った結果を得られることや、ULP の性質などに注意します。

まず、$\ulp(\sqrt{a}+\sqrt{b})\ge 4\ulp(\sqrt{c})$ の場合、 $$ \ulp(\sqrt{a}+\sqrt{b})\cdot 2^{p-1}\le \sqrt{a}+\sqrt{b} $$ より $$ \begin{aligned} \sqrtf(a)\oplus\sqrtf(b) &\ge \ulp(\sqrt{a}+\sqrt{b})\cdot 2^{p-1}-\tfrac12\ulp(\sqrt{a}+\sqrt{b}) \\ &\gt \ulp(\sqrt{a}+\sqrt{b})\cdot 2^{p-2} \\ &= \ulp(\sqrt{c})\cdot 2^{p} \\ &\ge \sqrtf(c) \end{aligned} $$ なので、$(\sqrtf(a)\oplus\sqrtf(b))-\sqrtf(c)\ge 0$ となります。

次に、$\ulp(\sqrt{a}+\sqrt{b})=2\ulp(\sqrt{c})$ の場合、 $$ \begin{aligned} \sqrtf(a)\oplus\sqrtf(b) &\ge \ulp(\sqrt{a}+\sqrt{b})\cdot 2^{p-1}-\tfrac12\ulp(\sqrt{a}+\sqrt{b}) \\ &= \ulp(\sqrt{c})\cdot 2^p - \ulp(\sqrt{c}) \\ &\ge \sqrtf(c) - \ulp(\sqrt{c}) \end{aligned} $$ なので、$(\sqrtf(a)\oplus\sqrtf(b))-\sqrtf(c)\ge -\ulp(\sqrt{c})$ となります。

最後に $\ulp(\sqrt{a}+\sqrt{b})=\ulp(\sqrt{c})$ の場合、

$$ \begin{aligned} \sqrtf(a)\oplus\sqrtf(b) &\ge \roundf(\sqrt{a}+\sqrt{b}) - \ulp(\sqrt{a}+\sqrt{b}) \\ &\ge \roundf(\sqrt{c}) - \ulp(\sqrt{a}+\sqrt{b}) \\ &= \sqrtf(c) - \ulp(\sqrt{c}) \end{aligned} $$ なので、こちらも $(\sqrtf(a)\oplus\sqrtf(b))-\sqrtf(c)\ge -\ulp(\sqrt{c})$ となります。

また、$\sqrtf(a)\oplus\sqrtf(b)\lt \sqrtf(c)$ の前提では $$\ulp(\sqrt{c})\cdot2^{p-1} \le \sqrtf(a)\oplus\sqrtf(b)\le \ulp(\sqrt{c})\cdot2^p$$ となります（刻み幅が $\ulp(\sqrt{c})$ の区間にあるということ）。

よって、差が $\ulp(\sqrt{c})$ 以下であることから、$(1-\tfrac12)\ulp(\sqrt{c})=\tfrac12\ulp(\sqrt{c})$ より真に大きい値を $\sqrtf(a)\oplus\sqrtf(b)$ に足すことで、$\sqrtf(c)$ 以上の値に丸められます。 真に大きくする理由は、丸めの tiebreak で意図しない方に丸められるのを防ぐためです。

$c = 10^9$ かつ $p = 64$ の下では、これを満たす最小の値は $2^{\floor{\log_2(\sqrt{10^9})}-64}\cdot(1+2^{-63})$ であり、 $$ {\footnotesize 8.881784197001252324} {\scriptsize 3520183169} {\scriptsize 2018042652} {\tiny 79889712924636592690508241076940976199693977832794189453125} \times 10^{-16} $$ です。

なお、リテラルとして書いた際にも正確に丸められるため、$2^{\floor{\log_2(\sqrt{10^9})}-64}\cdot(1+2^{-64})$ より少しでも大きい値であれば上記の値に丸められます。大きすぎると別の値に丸められるので注意です。

constexpr long double eps = 8.8817841970012523238705358308233714632639944856462318296345254120538470488099846988916397094726562500000000000000000000000001e-16L;

/* 短めに書きたい人向け */
constexpr long double eps = 8.881784197001252324e-16L;
constexpr long double eps = 8881784197001252324e-34L;

もちろん、何らかの事情がない限りは、ぎりぎりの値を設定する必要はないため、9e-16L などと簡単な値を使いたいです。 そのくらい大きくしても問題ないことを示すため、次の節が重要です。

$\varepsilon_U'$ の下界

$\sqrt{a}+\sqrt{b}\lt\sqrt{c}$ を考えます。 まず、誤差のない世界での計算を考えます。

$$ (\sqrt a+\sqrt b+\sqrt c)(\sqrt a+\sqrt b-\sqrt c)(\sqrt a-\sqrt b+\sqrt c)(\sqrt a-\sqrt b-\sqrt c) $$ は整数となります。各因子の符号は $$ \underbrace{(\sqrt a+\sqrt b+\sqrt c)}_{\gt0}\underbrace{(\sqrt a+\sqrt b-\sqrt c)}_{\lt0}\underbrace{(\sqrt a-\sqrt b+\sqrt c)}_{\gt 0}\underbrace{(\sqrt a-\sqrt b-\sqrt c)}_{\lt0} $$ であり、整数性と合わせて $1$ 以上であることがわかります。符号が負のものは絶対値を考えるとして、 $$ |\sqrt a+\sqrt b-\sqrt c| \ge \frac1{(\sqrt a+\sqrt b+\sqrt c)(\sqrt a-\sqrt b+\sqrt c)\cdot|\sqrt a-\sqrt b-\sqrt c|} $$ であり、右辺の分母を上から評価したいです。

まず、$(\sqrt a+\sqrt b+\sqrt c)\lt 2\sqrt c$ です。 次に、相加相乗平均の関係に気をつけつつ $$ \begin{aligned} (\sqrt a-\sqrt b+\sqrt c)\cdot|\sqrt a-\sqrt b-\sqrt c| &= -(\sqrt a-\sqrt b+\sqrt c)(\sqrt a-\sqrt b-\sqrt c) \\ &= -( (a+b-2\sqrt{ab})-c) \\ &= c - (a+b) + 2\sqrt{ab} \\ &\le c - (a+b) + (a+b) = c \end{aligned} $$ とわかります。よって、 $$ \sqrt c-(\sqrt a+\sqrt b) \gt \frac1{2c\sqrt c} $$ となります。すなわち、 $$ \sqrt a+\sqrt b\lt \sqrt c - \frac1{2c\sqrt c} $$ です。

$$ \begin{aligned} \sqrtf(a)\oplus\sqrtf(b) &\le \roundf(\sqrt{a}+\sqrt{b}) + \ulp(\sqrt{a}+\sqrt{b}) \\ &\le (\sqrt{a}+\sqrt{b})+\tfrac32\ulp(\sqrt{a}+\sqrt{b}) \\ &\lt \left(\sqrt c - \frac1{2c\sqrt c}\right) + \tfrac32\ulp(\sqrt{c}) \\ &= \sqrt c - \frac1{2c\sqrt c} + 3\cdot2^{\floor{\log_2(\sqrt{c})}-p}, \\ \sqrtf(c) &\ge \sqrt{c} - \tfrac12\ulp(\sqrt{c}) \\ &\ge \sqrt{c} - 2^{\floor{\log_2(\sqrt{c})}-p} \end{aligned} $$ なので、それらの差は $$ \begin{aligned} &\phantom{{}={}} \sqrtf(c) - (\sqrtf(a)\oplus\sqrtf(b)) \\ &\gt (\sqrt{c} - 2^{\floor{\log_2(\sqrt{c})}-p}) - \left(\sqrt c - \frac1{2c\sqrt c} + 3\cdot2^{\floor{\log_2(\sqrt{c})}-p}\right) \\ &= \frac1{2c\sqrt c} - 4\cdot2^{\floor{\log_2(\sqrt{c})}-p} \end{aligned} $$ となります。この差が $k\ulp(\sqrt{c})$ のとき、丸め誤差によって $$ (\sqrtf(a)+\sqrtf(b))\oplus\underbrace{(k-\tfrac12)\ulp(\sqrt{c})}_{\varepsilon} = \sqrtf(c) $$ になりえます。$\varepsilon$ を $(k-\tfrac12)\ulp(\sqrt{c})$ よりわずかに小さく設定すれば大丈夫そうです。

ということで、$c = 10^9$, $p = 64$ で計算してみると $$ \begin{aligned} \frac1{2c\sqrt c} - 4\cdot2^{\floor{\log_2(\sqrt{c})}-p} &= \tfrac12(10^9)^{-3/2} - 4\cdot 2^{\floor{\log_2(\sqrt{10^9})}-64} \\ &= \tfrac12\cdot 10^{-27/2} - 4\cdot 2^{-50} \\ &\approx 1.403 \times 10^{-14} \end{aligned} $$ となります。 $$ \underbrace{7.5\ulp(\sqrt{10^9})}_{\gt 1.332\times10^{-14}} \lt 1.403\times 10^{-14} \lt \underbrace{8\ulp(\sqrt{10^9})}_{\lt 1.422\times10^{-14}} $$ なので、$2k$ の整数性より $8\ulp(\sqrt{10^9})$ の差で下から抑えられそうです*4。 $7.5\ulp(\sqrt{10^9})$ 未満の最大の正規化数は $(1.875 - 2^{-63})\times 2^{16-63}$ であり、 $$ {\footnotesize 1.332267629550187848} {\scriptsize 43132080393349494087} {\tiny 7760882296602907258475934071384472190402448177337646484375} % {\footnotesize 1.332267629550187848} % {\scriptsize 4698394028} % {\scriptsize 2123965793} % {\tiny 88804411483014536292379670356922360952012240886688232421875} \times 10^{-14} $$ です。

こちらも、リテラルとして書いた際の丸めを考慮すると、$(1.875 - 2^{-64})\times 2^{16-63}$ より少しでも小さい値であれば上記の値に丸められます。

constexpr long double eps = 1.332267629550187848469839402821239657938880441148301453629237967035692236095201224088668823242187499999999999999999999999999e-14L;

/* 短めに書きたい人向け */
constexpr long double eps = 1.3322676295501878484e-14L;
constexpr long double eps = 13322676295501878484e-33L;

hack について

冒頭で、$8.5 \ulp(10^9)$ よりわずかに小さい値 1.5099033134902128948e-14L を使った解法は hack 可能と書きました。

hack できるケースの例は下記です。

239488342 239488343 957953370
247729794 247784320 991028225
248961062 248992620 995907363
249514080 249568802 998165761
249999995 249999996 999999982

これ以外にもたくさんあって、見つけた限りでも 649,152 個ありました。

ケースに関して

ローカル環境が M2 Mac なのですが、long double の精度が double と同じ $p=53$ でつらいです。 数分かかるコードを AtCoder のコードテストに投げるわけにもいかないので、Docker でそういう環境を作って対処しました。$p=113$ の環境も作ることができたので満足です。

$\sqrt a+\sqrt b=\sqrt c$ かつ $\sqrtf(a)\oplus\sqrtf(b)\lt\sqrtf(c)$ の下で、誤差が大きくなるケースについても探索しました。 $\sqrt{c}\in\N$ となるケースでは（制約の範囲では）誤差が出ないので、square-free な整数 $m\ge 2$ に対して $i_a\sqrt m+i_b\sqrt m=(i_a+i_b)\sqrt m$ のケースを全探索しました。

8 268609842 268702562
8 799840008 800000000
1002528 236661768 268470792
99904500 467737920 999980820
100017228 400692747 901090683
100026368 400218632 900407048
107717610 430870440 969458490
109996887 439987548 989971983
249961152 250015923 999954147

など、235,172,054 個のケースで、$2^{-49} \approx 1.776\times 10^{-15}$ の誤差を出すことができました。

上記の $\sqrt{a}+\sqrt{b}=\sqrt{c}$ となるケース (2,809,012,244 個) と、後述の $\sqrt{a}+\sqrt{b}\lt\sqrt{c}$ の差が極小となるケース (250,056,424 個) すべてについて、今回設定した eps が正しく判定できていることを確認しました。前者は 11 分、後者は 1 分程度で済みました。

差の下界について

$\sqrt{a}+\sqrt{b}\lt\sqrt{c}$ のケースにおける $\sqrt{c}-(\sqrt{a}+\sqrt{b})$ の下界についてです。

この節の内容は eps の設定には直接関係ない（範囲を絞るだけなら前の節だけで十分）ですが、自分で解決できなかったので載せておきます。数学に自信のあるフォロヮをお待ちしております。

$$ f(c) = \min_{\substack{\sqrt{a}+\sqrt{b}\lt\sqrt{c'}\\c'\le c\\a, b, c\in\N}} \sqrt{c}-(\sqrt{a}+\sqrt{b}) $$ を考えます。

Conjecture 1: $c\le 10^5$ くらいで実験したところ、$c$ を昇順に見ていったときに $f(c)$ の値が更新されるのは、$c$ が下記の形で書けるとき (if and only if) になりそうです。

$c$	$a$	$b$	$c\bmod 4$	$f(c)$
$2(2k+1)$	$k$	$k+1$	$2$	$\varepsilon_0(c)$
$(2k+1)(2k+3)$	$k(k+1)$	$(k+1)(k+2)$	$3$	$\varepsilon_1(c)$
$2(k+1)(k+2)$	$\tfrac12 k(k+1)$	$\tfrac12(k+2)(k+3)$	$0$	$\varepsilon_2(c)$
$(2k+1)(6k+3\pm2)$	$k(3k\pm1)$	$(k+1)(3(k+1)\pm1)$	$1$	$\varepsilon_3(c)$

ただし、$\varepsilon_i(x)$ は $x\ge i$ で定義される関数 $$ \varepsilon_i(x) = \sqrt{x} - \frac12\left(\sqrt{x+i-2\sqrt{ix+1}} + \sqrt{x+i+2\sqrt{ix+1}}\right) $$ です。

Conjecture 2: 任意の $i$ について $\varepsilon_i(x) \sim \tfrac12x^{-3/2}$ が成り立ちそうです。

Wolfram|Alpha 的には合っていそうでしたが、数学力が足りず、自分では示せませんでした。

Conjecture 3: また、任意の $i$ について $x\ge i+1$ ならば $\varepsilon_i(x) \lt \varepsilon_{i+1}(x)$ が成り立ちそうです。

グラフをプロットした感じは合っていそうでしたが、手計算はうまくいかず、投げ出してしまいました。

まとめ

Sqrt Inequality に関して、$c\le 10^9$ かつ $p=64$ であれば $$(1+2^{-63})\times 2^{-50} \le \text{\texttt{eps}} \le (1.875 - 2^{-63})\times 2^{-47}$$ の範囲で設定すればよさそうということがわかりました。10 進 4 桁だと $$8.882\times 10^{-16} \le \text{\texttt{eps}} \le 1.332\times 10^{-14}$$ となります。

また、ジャッジ上で AC となるコードを hack できうるケースをたくさん見つけることができました。

浮動小数点数の誤差の評価が大変ということもわかりました（実際にはもっとゆるっと評価していい部分をがんばって詰めようとしたため、必要以上に大変になったという部分はあります）。これからもやっていきたいと思いました。

参考

twitter.com

こういうの、下からの評価の方が上からの評価よりも難しくなりやすく。逆数を上から評価しましょう。

x_0 = \sqrt{a}+\sqrt{b}+\sqrt{c}
x_1 = \sqrt{a}+\sqrt{b}-\sqrt{c}
x_2 = \sqrt{a}-\sqrt{b}+\sqrt{c}
x_3 = \sqrt{a}-\sqrt{b}-\sqrt{c}

とすると、x_0x_1x_2x_3 は整数になります。

— maspy (@maspy_stars) 2020å¹´3月14æ—¥

このツイートでは $$\sqrt{c}-(\sqrt{a}+\sqrt{b})\ge \tfrac1{27}\cdot (10^9)^{-3/2}$$ だけを示していますが、方針で詰まっているときに読んでとても参考になりました。やっぱり数学に明るくて経験もある方は違うなと思いました。

• ABC 205 A
• ABC 263 Ex
• ABC 286 Ex

ぽかいこちゃんリスペクト。

感情

浮動小数点数の誤差については未証明でもおっけーみたいな風潮ありますよね。 未証明貪欲についてあれこれ言っている人が double で適当に計算して AC することに対して平気なのはダブルスタンダードなのではないでしょうか。

実用的にやばい例をまだ思いつけていないのですが、ある述語 $P(x)$ が、無限精度で計算したときは「$x\in(-\infty\ldots \alpha]$ では true となり、$x\in(\alpha\ldots{+\infty})$ では false となる」という条件を満たすとしても、$P(x)$ を浮動小数点数型で計算した場合にその性質が満たされるとは限りません。

たとえば $x \le 0 \vee (\tfrac 1x)\cdot x = 1$ という条件は上記の性質を満たします ($\alpha=+\infty$) が、浮動小数点数型では多くの反例があります。

$p$	型	反例
$24$	float (32-bit)	$41$
$53$	double (64-bit)	$49$
$64$	long double (80-bit)	$41$
$113$	long double (128-bit)	$43$

また、$[1\ldots 2)$ の範囲の最小の反例は、上記の $p$ の順に次のようになります。$p=113$ については探索中です*5。 $$ \begin{aligned} \texttt{801467}_{(16)}\times2^{-23} &= {\footnotesize 1.000622630119323730}{\scriptsize 46875}, \\ \texttt{1000000f5cbf2a}_{(16)}\times2^{-52} &= {\footnotesize 1.000000057228997096}{\scriptsize 81428248586598783731}{\tiny 4605712890625}, \\ \texttt{800000005a82799a}_{(16)}\times 2^{-63} &= {\footnotesize 1.000000000164636126}{\scriptsize 99697816044164255799}{\tiny 842067062854766845703125}. \end{aligned} $$

$p = 113$ は、全探索だと無理そうだったので、考えて求めました。別途記事にします。$1+(2^{57}+\varepsilon)\cdot 2^{-112}$ くらいだったので、$1$ から昇順に求めるのをやめて正解でした。 $$ \begin{aligned} &\phantom{{}={}} \texttt{1000000000000022df0668c89bf13}_{(16)}\times 2^{-112} \\ &= {\footnotesize 1.000000000000000030}{\scriptsize 24593730708203819677}{\tiny 89344748946071241091121526548576685378133532822175766341388225555419921875}. \end{aligned} $$

C/C++ で計算すると、float が勝手に double に変換される (floating-point promotion) ので誤差が生じないように見えましたが、$p=24$ で計算していないので不正ですね。 なお、各値が正規化数となる範囲においては、上記のように $(1.0\oslash x)\otimes x = 1-2^{-p}$ になることはあるものの、上方向にズレて $1+2^{1-p}$ とはならないことを示せます（別の記事に書きそう）（非正規化数についてはまだ考えていません）。

↓ 追記：上記二つのトピックに対して書きました ↓

rsk0315.hatenablog.com

二分探索で何かしらの計算をしていて、的外れな値で意図しない true/false が返ってきて変なことになる可能性は全然あると思っていて、そうしたことに無頓着なのはどうなの？という気が最近しています。 こわいので、最近は、可能な限り有理数で二分探索するようにしています。

atcoder.jp

とはいえ、計算途中でオーバーフローしないことの担保すら “非本質な細部” のように思っている人も多数いそうな感じがしますし、要求しすぎでしょうか。 実際、これらの評価をコンテスト中にやれと言われたら、それが楽しいとはあまり思わなさそうです（というか時間が足りなさそう）（別のジャンルの競技としてはありかもしれませんが）。

それはそれとして、元の問題の解説で「整数でやりましょう」で済ませるのではなく「$\varepsilon=10^{-14}$ だとうまくいくことを説明できる」と書いてあるのがやっぱりすごいなと思いました。

ABC 169 C なども誤差に関する教育的な問題です。コンテスト当時は long double で解くだけで AC になってしまうようでしたが、現在はケースが追加されているようです。こちらについても誤差評価を考えてみてはいかがでしょうか。

おまけ

Python の Decimal で精度を上げて計算した場合も、同様にたくさんのケースで hack が可能です。

atcoder.jp

2 2 8
8 8 32
2 8 18

この子たちは、小さいケースながら、多くの精度・丸め設定において誤判定させることができました。

おわり

お疲れさまでした。

*1:複数ある場合は仮数部が偶数であるものとします。いわゆる roundTiesToEven です。

*2:$\roundf(x+y)$ と $\hat x\oplus\hat y$ の ULP は異なるケースがありますが、ここでは小さい方を採用するとします。

*3:多くの問題は、このような “細部” に気を払わずに解けるようになっていますが、そのせいで「これは職人がノリで決める値だ」「よくわからないが適当に書いておけばよい」と思考停止している人が多いのはあまり好ましいことではない気がします。

*4:$2\ulp(\sqrt{a}+\sqrt{b})=\ulp(\sqrt{c})$ となるケースで $k$ は半整数となることがありそうです。

*5:$1.9\times 10^{12}$ 個調べましたがまだ出ていません。見つかり次第追記します。もしかしてこれ $1+2^{-56}$ くらい未満には存在しなかったりしますか？