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

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

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

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 &= (734362212, 578450837, 864377693, 680862920), \\ (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 7.69\times10^{-9}$	$\approx 4.92\times 10^{-7}$	$\approx 4.85\times 10^{-7}$	$\approx63$
$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
734362212 578450837
864377693 680862920

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:未来のえびちゃん見てますか 👋