N - 有限べき級数/Finite Power Series Editorial by rsk0315

誤差で WA になる解法とそのケース

公式解説 の記号を流用します。

考察

\[ (f_l\circ f_{l+1}\circ\dots\circ f_r)(0) = (f_l\circ f_{l+1}\circ\dots\circ f_N)(0) - (f_{r+1}\circ f_{r+2}\circ\dots\circ f_N)(0)\cdot x^{r-l+1} \]

が成り立ちます。

なので、\(y_i = (f_i\circ f_{i+1}\circ\dots\circ f_N)(0)\) を求めておけば、\(y_l - y_{r+1}\cdot x^{r-l+1}\) として求められそうです。

解法

64-bit 浮動小数点型の加減算と乗算(誤差を含む)をそれぞれ \(\oplus\), \(\ominus\), \(\otimes\) と書くことにします。 また、\(x\), \(x^k\), \(y\) を 64-bit 浮動小数点型で表そうとした値をそれぞれ \(\hat x\), \(\hat p_k\), \(\hat y\) と書くことにします。 ここで、\(\hat p_k\), \(\hat y\) に関しては、真の値 \(x^k\), \(y\) を正確に丸めたものではなく、下記のアルゴリズムによって求めたものとします。\(\hat x\) に関しては、入力で受け取った \(x\) を正確に丸めたものとします。

  • 前計算:
    • \(\hat y_{N+1} = 0\) および \(\hat y_i = (\hat y_{i+1}\otimes \hat x) \oplus A_i\) として \(i = N, N-1, \dots, 1\) の順に求める。
    • \(\hat p_0 = 1\) および \(\hat p_{i+1} = \hat p_i \otimes \hat x\) として \(i = 1, 2, \dots, N\) の順に求める。
  • クエリ処理:\(\hat y_l\ominus (\hat y_{r+1}\otimes\hat p_{r-l+1})\) を出力する。

note: 浮動小数点型の四則演算は結合法則が成り立たないことなどから、アルゴリズムの細部によって計算結果が異なる場合があります。同様の考察に基づく解法でも、括弧のつけ方が一つ違うだけで出力が変わることは普通に起こります。

提出 #62630164

撃墜

さて、この解法で WA となるようなケースを構築します。

一般に浮動小数点型の減算 \(a\ominus b\) において、\(a\)\(b\) が近い値であり、かつ \(a\)\(b\) が誤差を含む値であるとき、計算結果の誤差が非常に大きくなり得ます。(桁落ちや catastrophic cancellation などと呼ばれます。)

そこで、\(\hat y_l\)\(\hat y_{r+1}\otimes \hat p_{r-l+1}\) がそれらの条件を満たすようなケースを考えます。

note: \(A\) の形を考えた後は、\(N\)\(x\) を全探索して十分な誤差が出るケースを探しています。

ケース 1

  • \(N = 36\)
  • \(A = (0, 1, 10^9, 10^9, \dots, 10^9)\)
  • \(x = 0.9\)
  • \(Q = 1\)
  • \((l_1, r_1) = (1, 2)\)

このとき、真の値は \(0+1\cdot x = 0.9\) ですが、上記のアルゴリズムでは \(0.9-1.335\cdot 10^{-6}\) 程度の値を出力します。真の値などは下記の通りです。

\[ \begin{aligned} y_1 &= x + 10^9\cdot(x^2+\dots +x^{N-3}) \\ &= x + 10^9\cdot x^2\cdot\frac{1-x^{N-2}}{1-x} \\ &= 7874716005.{\small 450608255881598}{\footnotesize 52125227359}, \\ y_3\cdot x^2 &= 10^9\cdot x^2\cdot\frac{1-x^{N-2}}{1-x} \\ &= 7874716004.{\small 550608255881598}{\footnotesize 52125227359}, \\ \hat y_1 &= 7874716005.{\small 450611114501953}{\footnotesize 125}, \\ \hat y_3\otimes \hat p_2 &= 7874716004.{\small 550612449645996}{\footnotesize 09375}, \\ \hat y_1 \ominus (\hat y_3\otimes \hat p_2) &= 0.{\small 899998664855957}{\footnotesize 03125}. \end{aligned} \]

制約の範囲内で、\(x\lt 1\) のとき \(y_3\cdot x^2 \approx 1.999\cdot 10^{14}\) 程度には大きくできるため、catastrophic cancellation によって \(1.999\cdot 10^{14}\times 2^{-52} \approx 4.439\cdot 10^{-2}\) 程度の絶対誤差が生じうることは予想できます。実際、\((N, x) = (199983, 0.999999)\) のとき \(3.125\cdot 10^{-2}\) 程度の誤差が生じました。

ケース 2

  • \(N = 55\)
  • \(A = (0, 0, \dots, 0, 10^9)\)
  • \(x = 0.999999999 = 1 - 10^{-9}\)
  • \(Q = 1\)
  • \((l_1, r_1) = (1, N-1)\)

このとき、真の値は \(0\) ですが、上記のアルゴリズムでは \(1.073\cdot 10^{-6}\) 程度の値を出力します。

\[ \begin{aligned} y_1 &= y_N\cdot x^{N-1} \\ &= 10^9\cdot x^{N-1} \\ &= 999999946.{\small 000001430999975}{\footnotesize 196000316250996}{\scriptsize 837490025827164}{\tiny 822899441040465}{\tiny 78468206376}, \\ \hat y_1 &= 999999946.{\small 000002622604370}{\footnotesize 1171875}, \\ \hat y_N\cdot \hat p_{N-1} &= 999999946.{\small 000001549720764}{\footnotesize 16015625}, \\ \hat y_1 \ominus (\hat y_N\cdot \hat p_{N-1}) &= 0.{\small 000001072883605}{\footnotesize 95703125}. \end{aligned} \]

note: \(\hat p_i\) を漸化式ベースで求めるのではなく pow などの関数を用いた場合でも、(pow の実装によりますが)\((N, x) = (51, 0.9999)\) などのケースで許容できない誤差が生じます。

本質的には、\(a_N\otimes (x\otimes\dots\otimes x)\)\((a_N\otimes x)\otimes\dots\otimes x\) が等しいとは限らないことに起因します。\(k\) 回の乗算では \((1+2^{-53})^k-1\) 程度の相対誤差が生じうるため、制約の範囲内では \((1+2^{-53})^{199999}-1 \approx 2.220\cdot 10^{-11}\) 程度になります。\(x\lt 1\) の範囲では \(a_N\cdot x^{N-1} \approx 9.999\cdot 10^8\) なので、これらを掛けて \(2.220\cdot 10^{-2}\) 程度の絶対誤差が生じうることは予想できます。

\((N, x) = (199999, 1-10^{-15})\) のケースでは、\(9.105\cdot 10^{-3}\) 程度の絶対誤差が生じました。\(x\) の桁数が少なめのケースとしては、\((N, x) = (294, 0.999)\) のとき \(1.192\cdot 10^{-6}\) 程度の絶対誤差が生じました。

コピペ用

36
0 1 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000 1000000000
0.9
1
1 2
55
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1000000000
0.999999999
1
1 54

あとがき

浮動小数点型で減算をするような解法では、今回のように catastrophic cancellation させられるような入力が存在することが多いです。 とはいえ、個別の解法たちに対して撃墜ケースを用意するのは大変であることも多いため、そうした解法でも AC になってしまうことはしばしばありそうな気はします。

コンテスト後には、「本当に正しい解法なのか?」ということを確認してみるのも重要だと思います。

posted:
last update: