不必要に浮動小数点数を使うの怖い

競プロにおいて、「何かしらの整数を大小比較して、左辺が大きければどうのこうの」みたいなことをしたい局面はそれなりにあります。 特に、境界付近でも正確に評価できる必要がある状況も多いです。 そうした状況で浮動小数点数を使うのはできるだけ避けたいですという気持ちがあります。

å°Žå…¥

まず、浮動小数点数は誤差が出ます*1。

初心者のうちは「浮動小数点数では \(1\div 10\) を正確には計算できない」と聞いても、「たとえば

assert_eq!(1.0 / 10.0, 0.1);

とかはエラーにならないし、問題なく計算できているのでは？」みたいな誤解をすることもあります。 実際には、0.1 と書いた部分も計算前の時点で 0.1000000000000000055511151231257827021181583404541015625 のような値になっており、所望の状態にはなっていません。 1.0 / 10.0 もそういう値になっていたためエラーにならなかっただけですね。

• そもそも表せない実数が（たくさん）ある
• 表せる数同士の演算結果でも表せない数になることも（よく）ある

というのをちゃんと意識する必要があるかもしれません。 もちろん、表せなかった数同士で計算した結果がたまたま無限精度での計算結果と同じになることもあります。

assert_eq!(0.1 * 10.0, 1.0);

実際には、0.1000000000000000055511151231257827021181583404541015625 * 10.0 を計算しようとして、\(1.000000000000000055511151231257827021181583404541015625\) を表せずに 1.0 に丸められた結果、見かけ上は等しくなっているわけですね。

というわけで 1.0 / 10.0 * 10.0 == 1.0 は true ですが、たとえば 1.0 / 49.0 * 49.0 == 1.0 は false になります。 左辺は 0.99999999999999988897769753748434595763683319091796875 のような値になっています。

また、他によくある話として 0.1 + 0.2 != 0.3 というのもあります。 誤差をどう丸めるかにはいくつかのモードがありますが、「表せる数のうち最も近いものに丸める。ちょうど中間にある場合は、最も近い偶数の方向に丸める」という（よくある）ものを仮定すると、

0.1 => 0.1000000000000000055511151231257827021181583404541015625
0.2 => 0.200000000000000011102230246251565404236316680908203125
0.3 => 0.299999999999999988897769753748434595763683319091796875

となりますが、加減乗除については「無限精度での演算をしてから丸めた値と等しくなる」というように規定されているので、0.1 + 0.2、すなわち 0.100000000000000005551... + 0.200000000000000011102... を無限精度で計算してから丸めた 0.3000000000000000444089... のようになり、0.3 を丸めた 0.299999999999999988897... とは等しくなりません。「0.1 + 0.2 と書いたんだから、2 進法とかにはせずに書いた時点の値を無限精度で計算してから丸めてくれ」というのはきつそうです*2。

0.1 + 0.2 => 0.3000000000000000444089209850062616169452667236328125

型について

double ã‚„ f64 と言われる型でどのように数を表すかについて触れます。 IEEE 754 という規格で binary64 として定義されているものです*3。 実数の他に \(\infty\) に対応する値 Infinity や、非数（定義域外の値を関数に入力したときの返り値になったりする）に対応する値 NaN などもこれらの型には含まれますが、一旦それらは無視します*4。

例として 3.141592653589793115997963468544185161590576171875 を挙げます。これは以下のようなビット列で表されます。

0 10000000000 1001001000011111101101010100010001000010110100011000

三つに分けて書いた部分をそれぞれ \( (s, e, b_{51} b_{50} \dots b_1 b_0)\) とおくと、次の式で表される値に対応します。 \[ (-1)^s \cdot 2^{e-1023} \cdot (1.b_{51} b_{50} \dots b_1 b_0)_{(2)}. \] つまり、以下のような感じです。 \[ \begin{aligned} {} &{\phantom{{}={}}} 2^{1024-1023} \cdot 1.1001001000011111101101010100010001000010110100011000_{(2)} \\ &= 2\cdot 1.1001001000011111101101010100010001000010110100011000_{(2)} \\ &= 2\cdot (1 + 2^{-1} + 2^{-4} + 2^{-7} + 2^{-12} + \dots + 2^{-48} + 2^{-49}) \\ &= 3.141592653589793115997963468544185161590576171875 \end{aligned} \]

というわけで、一番上の \(1\) がある桁から見て \(52\) 桁以上離れた \(1\) があるような数 (\(2^k\cdot 1.\underbrace{00\dots 0}_{52}1_{(2)}\)) や、十分大きい数 (\(2^{1024}\)) などは表せないことがわかります。このことから当然、有限小数でない実数も表せないことがわかります。\(2\) 進法において \(0.1_{(10)}\) は有限小数ではないため、\(0.1\) も表せません。

bartaz.github.io

こういうサイトがあって、オサレな UI で遊べます。

有限小数について

\(b\) 進法で既約分数 \(p/q\) が有限小数となる条件について考えておきます。簡単のため \(0\lt p/q\lt 1\) とします。

ある有限の \(k\) について \[ p/q = (0.x_1 x_2\dots x_k)_{(b)} \] とすると、両辺を \(b^k\) 倍して \[ b^k\cdot p/q = (x_1x_2\dots x_k)_{(b)} \] です。すなわち、\(b^k\cdot p/q\) が整数になる必要があります。\(p\) と \(q\) は互いに素なので \(b^k/q\) が整数になる必要があります。 このとき、（\(b\) と \(q\) の各素因数について約分していくみたいなことを考えると）\(q\) の素因数を \(b\) がすべて持っている必要があります。 逆に、\(q\) の素因数を \(b\) がすべて持っていれば十分です。

\(0.1 = 1/10\) に関しては、\(10\) の素因数 \(5\) を \(2\) が持っていないため、\(2\) 進法では表せないわけです。

嫌な例

単純な分数の比較

「\(x / y \le z / w\) か？ (\(x, y, z, w\in \Z\cap[-10^2, 10^2]\), \(0\notin\{y, w\}\))」というのを知りたいケースはまぁよくあります。 これは、分母を払って \(xw\le yz\) とすれば整数で比較できます*5が、浮動小数点数で計算してしまう初心者は多いです。

ざーーっと調べてみた感じ、hack できるケースは見つからなかったのですが、コンテスト本番で不安になりながら書きたくはないので、整数でやりたいです。

大きい整数の比較

「\(xy \ge z\) か？ (\(x, y, z\in\Z\cap[0, 10^{18}]\))」というのを知りたいケースはそこそこあるでしょう。x * y がオーバーフローするので整数でやりたくないと思ってしまう人が多そうです。

fn f(x: i64, y: i64, z: i64) -> bool {
    let x = x as f64;
    let y = y as f64;
    let z = z as f64;
    x * y >= z
}

たとえば、\(10^{18}-1 = 27027027\times 37000000037\not\ge 10^{18}\) ですが、f(27_027_027, 37_000_000_037, 10_i64.pow(18)) は true になってしまいます。 左辺が誤差で 1_000_000_000_000_000_000.0 になってしまい、等しくなってしまうんですね。

参考袋 を見つつ、\(x \ge \ceil{z/y}\) と変形します。\(y = 0\) に気をつけて、

fn f(x: i64, y: i64, z: i64) -> bool {
    if y == 0 { z == 0 } else { x >= (z + y - 1) / y }
}

のようになる気がします。

もっとも、i128 のような大きい整数型を使うという手もありますが、制約が \(10^{36}\) とかだったら？と考えると趣旨は伝わると思います。

小数点の移動

「小数点以下高々 \(9\) 桁の実数 \(0\lt x\lt 10^4\) が与えられます。\(10^9\cdot x\) を出力してね」というの を考えます。

適当に 1.0e9 倍して少し足しておけばいいでしょうみたいな気持ちになる人も多いと思うのですが、

fn f(x: f64) -> i64 { (1.0e9 * x + 1.0e-6).floor() as i64 }

などと書くとうまくいきません。"32.000000001" が与えられたとき、これをそのまま f64 に変換すると 32.00000000099999653002669219858944416046142578125 となり、32000000000.999996185302734375 となり、32000000000 となります。

文字列としてがんばって処理するか、.round() as i64 のようにするのが無難そうです*6。

平方根の比較

これは有名。「\(\sqrt{a}+\sqrt{b}\lt\sqrt{c}\) ですか？」というやつです。

解説 を読んだり、maspy さんのツイート を読んだり。

累乗関連

Rust だと整数にも .pow() が用意されているのであんまりなさそうですが、C++ だと浮動小数点数にしか用意されていないせいか、整数の累乗をしたいときに浮動小数点数を介してしまう人が多いです。

「集合 \(S, T\) が与えられます。\(|S\cap T|\) は？ (\(S, T\subseteq\{0, 1, \dots, 62\}\))」というのがあったとき、

// C++
long long x = 0, y = 0;
for (auto si: s) x += pow(2, si);
for (auto ti: t) y += pow(2, ti);
return __builtin_popcountll(x & y);

みたいなことを書いてしまう人もいます。pow が double を返すので、以下のような変換が暗黙になされています。

// x += pow(2, si);
x = (long long)((double)x + pow(2, si));

その結果、たとえば \(S = \{0, 53\}\) のようなとき、期待する値が \(9007199254740993\) なのに対し、誤差の関係で x == 9007199254740992 になってしまいます。 x |= 1LL << si のように書くのが無難な気がします。Rust においては x |= 1 << si と書いても問題ないです。

Fibonacci æ•°

\(n\) 番目 (\(1\le n\le 10^{18}\)) の Fibonacci 数を \(998244353\) で割ったあまりを求めてね、と言われて、\(F_n = (\sqrt{5})^{-1}\cdot (\phi^n-(-\phi^{-1})^n)\) を浮動小数点数演算で求めようとする初心者もいる気がします。

たとえば \(F_{72} = 498454011879264\) ですが、上式をそのまま計算すると、498454011879265.2 とかになってしまいそうです。

こっち はギャグの記事。

もっと単純な比較で hack できそうな例を出して初心者を脅かしたかったんですが、いいのが思いつきませんでしたね。

sqrt

\(0\le n\le 10^{18}\) で \(\floor{\sqrt{n}}\) を求めてね、と言われて、std::sqrt(n) と書くとだめです。 sqrt が悪いというよりは、double を受け取る関数に long を渡した際に暗黙に型変換が起きているのが原因なのですが。

たとえば \(n = 2^{54}-1\) のとき、double では \(2^{54}\) に丸められ、\(2^{27}\) が返されたりします。実際には \(\floor{\sqrt{n}} = 2^{27}-1\) なので嘘になりますね。 平方数より少し小さい整数で、上に丸められてその平方数以上になるケースであれば、なんでも hack できると思います。

指数探索なり二分探索なり Newton 法なり整数でやる方針にするか、double で返ってきた値の周辺を線形探索する方針がよさそうです。

おまけ

これは「整数を使えばいいですよ」系とは違って、「こういう例もあります。嫌ですね」みたいなのを紹介するコーナーです。

誤差のせいで定義域を超えるやつ

敢えて作ったような例ですが、「\(x\) が与えられます。\(\sqrt{\sin(\pi x)}\) を出力してね」という問題を考えます。\(x = 2\) のときなどにうまくいきません。

use std::f64::consts::PI;

fn main() {
    println!("{:.50}", (2.0 * PI));
    println!("{}", (2.0 * PI).sin());
    println!("{}", (2.0 * PI).sin().sqrt());
}

以下のような出力になります。

6.28318530717958623199592693708837032318115234375000
-0.00000000000000024492935982947064
NaN

0.0 になってくれるでしょというのは思い込みなわけですね。

\(2\pi = \tau\) を表す定数 TAU も用意されていますが、TAU.sin().sqrt() ã‚‚ NaN になります。

他にも、asin (\(\arcsin\)) の引数が 1.0 よりわずかに大きかったときも NaN になったりします。 定義域が \( (-\infty, \infty)\) でない関数*7を使うときは特に気をつけましょう。

意図しない感じに丸められるやつ

これも敢えて作った例ですが「\(1\le x\le 9007199254740995\) が与えられます。\(1/(9007199254740996-x)\) を出力してね」という問題を考えます。 x = 9007199254740995.0 とすると、x == 9007199254740996.0 となるように丸められ、1.0 / 0.0 == Infinity となります。

細かいやつ

\(1 / 3\) は 2 進法で有限小数にはならないので、\(\sqrt[3]{x} = x^{1/3}\) と思って計算すると少しずれることがあります。

println!("{:.050}", 5.0_f64.cbrt());
// 1.70997594667669705614798658643849194049835205078125

println!("{:.050}", 5.0_f64.powf(1.0 / 3.0));
// 1.70997594667669683410338166140718385577201843261719

競プロの範囲でこれが問題になることは少ないかもですが、頭の片隅に置いておいてもよいでしょう。

丸めに関して

加減乗除については、無限精度の演算結果を丸めたかのように振る舞ってくれます。

FMA (fused multiply-add; \(xy+z\)) というのがあって、丸めを一度しか行わないものとして定義されていますが、正しくない結果を返す実装も あるらしい です。

それ以外の関数 (\(\sin(x)\), \(\sqrt{x}\), \(x/y\), etc.) についてはどうでしょうか。 IEEE 754-2019 では required operations と recommended operations という枠が分かれていて、加減乗除などは前者、\(\sin\) などは後者に含まれています*8。 後者も無限精度の演算結果を丸めたかのようにするべきと書かれていますが、required ではないのでゆるっとした感じがする気がします。実際、GCC は一部を除いて特にちゃんとすることは目指していない らしいです。

table maker’s dilemma とかを調べるとよいかも。Lindemann の定理 とかも関連していそう。 関連してそうな資料、これ とか これ とか。

詳しい人のツイート：

C言語の場合はsinとかの数学関数がそもそもIEEE準拠かという問題があって、IEEE754を見に行く前にCの規格を確認するべき、みたいな話がある（Cの処理系がAnnex F準拠ならsqrtとかremainderは正しく丸められることが保証されるけど、sinは保証されない） https://t.co/Woml3ZkagN

— mod_poppo (@mod_poppo) 2022å¹´6月4æ—¥

その他

もしかすると x < y + eps (eps = 1e-8) みたいにするといいとかを聞いたことがある人もいるかもしれません。 それはどうしても浮動小数点数を使う必要がある場面で適切に使うべきであって、整数に帰着できる問題ではやらない方が無難そうです。

根本的に間違った解法の場合でも「eps を変えれば通るかも」などと言って無駄に誤答を出すことにもつながりがちです。

また、幾何や幾何以外の問題などで分数を扱う場合においても、適当な素数 \(p\) を用いて、\(x/y\) を \(x\cdot y^{-1}\bmod p\) として表すとよいケースもあります。(問題 / 解説) (問題 / 解説)

こと浮動小数点数においては、初心者のうちに「なんかこれでうまくいきそう（いくつかのケースではうまくいった）」と思ったものには、反例がたくさんある気がします。 そうした方法を使って祈るよりは、（十分に勉強して祈らず使うか）整数に帰着する方法を考えるのが賢明そうです。

関連問題

浮動小数点数を避けてやりたいな〜となる問題を挙げておきます。浮動小数点数でやって問題があるかどうかは気にしていないので、問題がない場合でも整数でやる練習としてやるとよいでしょう。

• ABC 169 C
  • hack
• ABC 226 A
• ABC 238 A
  • 入力が 1e160 くらいまである場合とかを考えてみてください。
• ABC 191 D
• ABC 205 A
  • 解説

他にもいい感じのがあったら募集しています。

計算に関して

docs.rs

浮動小数点数の数学関数の実装が載っています。 ハードウェアのサポートを得られない状況を想定して書かれていそうで、多くの状況ではハードウェア命令で処理されていそうな気がします。

近似とかをがんばったりしていそうです。

この辺は「不必要に〜」の文脈とは違いそうですが、ついでに載せちゃいました。

便利ちゃん？

誤差が減るように式変形してくれるらしいです。

herbie.uwplse.org

おわり

おわりです。

*1:誤差以外もいろいろ出そうと思えば出るかもしれません。

*2:コンパイル時の最適化によっては起きうる？ オプションとかにもよるかも。C/C++ ではコンパイル時と実行時で精度や丸め方が変わるのを許容するみたいな話がどこかにあったような気がするのでまた調べます。最適化で誤差関連の挙動を変えてもいいよーみたいなオプションはあったはず。

*3:Rust はそうっぽいけど C++ は必ずしもそうではないかも。GCC ではサポートしていそう。そうじゃない処理系があるかは知らないけど規格上は許容されていそう。

*4:さらに、非正規化数というのもありますが、無視します。

*5:というのは嘘で、\(yw\) の符号に気をつけましょう。

*6:反例を見つけたら教えてください。

*7:binary64 は無限大たちを含むので \([-\infty, \infty]\) かも。

*8:古いのは確認していないのですが、当時はなかったらしい です？ たぶん。