floor(n / 10.0) と floor(n * 0.1) は等価？

これを double で計算したときに常に等しくなりますか？というクイズです。 n は double で表せる整数であることにします。

答え

>>> 6755399441055749 / 10
675539944105574.9
>>> 6755399441055749 * 0.1
675539944105575.0

より、前者の floor と後者の floor は異なることになります。 すなわち、等価ではありません。

解説など

真であるなら証明をする必要があるんですが、偽だったので反例を挙げて終了です。 それはそれとして、「これの範囲では成り立つ」という話をしても損はしないので、してみましょう。

浮動小数点数型についてある程度の知識は持っているものとします。そうでない人のための記事は近々（次か次の次くらいに）書く予定です。

$\gdef\fround#1{(\hspace{-.2em}[#1]\hspace{-.2em})}$

まず、0.1 と言っているのは、double（かつ、多くの処理系でデフォルトになっている丸め方向）の文脈においては $$ \begin{aligned} &\phantom{{}={}} {\footnotesize 0.100000000000000005}{\scriptsize 55111512312578270211}{\tiny 81583404541015625} \\ &= \texttt{ccccccccccccd}_{(16)} \times 2^{-55} \\ &= (\tfrac45(2^{52}-1)+1)\cdot 2^{-55} \\ % &= (\tfrac45(2^{-3}-2^{-55})+2^{-55}) \\ &= \tfrac1{10}+\tfrac15\cdot 2^{-55} \\ \end{aligned} $$ と等しい値になります*1。

実数 $x$ を浮動小数点数型（ここでは double）で表せる値に丸めたものを $\fround{x}$ と書くことにします*2。 たとえば、$\fround{0.1} = 0.1 + \tfrac15\cdot 2^{-55}$ ã‚„ $\fround{10} = 10$ です。

また、浮動小数点数 $x$, $y$ に対して $\fround{x\times y}$ および $\fround{x\div y}$ を、それぞれ $x\otimes y$ および $x\oslash y$ と書くことにします。 浮動小数点数型の値同士の演算のみしか考えていないことに注意してください。

さて、今回やりたいのは $\floor{n\oslash 10}$ と $\floor{n\otimes\fround{0.1}}$ の比較です。

note: double の仮数部は $p = 53$ bits なので、$2^{53}$ 以上の（浮動小数点数型の）値 $x$ は整数となるため、$\floor{x}=x$ となります。

実数 $x\ge 0$ と整数 $k$ に対して $\floor{\fround{x}} = k$ となる条件は、$\fround{x}$ の整数部分が $k$ となることです。丸め誤差の関係で、$x$ の整数部分が $k$ であることとは同値ではありません。

さて、浮動小数点数で表せる整数 $n = 10q+r$ ($q\in\N, 0\le r\le9$) を考えます。$2^e\le 10q+r \lt 2^{e+1}$ とします。$0\le e\lt p$ としておきます。

$\fround{0.1} = \tfrac1{10} + \tfrac15\cdot 2^{-55} = \tfrac1{10}\cdot(1+2^{-54})$ なので、 $$ \begin{aligned} (10q+r)\cdot\fround{0.1} &= (q+\tfrac r{10})\cdot(1+2^{-54}) \\ &= q+(\tfrac r{10} + (q+\tfrac r{10})\cdot 2^{-54}) \end{aligned} $$ となります。一旦 $r$ が最大の $9$ である場合について考えます。 $$ \begin{aligned} (10q+r)\cdot\fround{0.1} &\le q+(\tfrac9{10} + \tfrac1{10}(10q+r)\cdot 2^{-54}) \\ &\lt q+\tfrac9{10} + \tfrac1{10}\cdot 2^{e+1-54} \\ &= q+\tfrac9{10} + \tfrac1{10}\cdot 2^{e-53}. \end{aligned} $$ $\tfrac1{10}(10q+r)\lt\tfrac18(10q+r)\lt 2^{e+1-3}$ より、$\tfrac1{10}(10q+r)$ の 0.5 ULP は $2^{e-56}$ 以下です。 小数部分を考えて、$\tfrac9{10}+\tfrac1{10}\cdot 2^{e-53} \lt 1-2^{e-56}$ の範囲においては、$(10q+r)\cdot\fround{0.1}$ は $q$ の方に丸められます。 これにより、$e\le51$ の範囲では $\floor{n\oslash 10} = \floor{n\otimes\fround{0.1}}$ となります。 また、同様の考察により、$r\le 8$ では $e=52$ でも $\floor{n\oslash 10} = \floor{n\otimes\fround{0.1}}$ となることが示せます。

よって、$\floor{n\oslash 10} \ne \floor{n\otimes\fround{0.1}}$ となる最小の $n$ は、$e = 52$ かつ $r = 9$ においてこれを満たす最小の $q$ によって作られることがわかります（この時点では存在性は示していませんが、これから作るので大丈夫です）。 $(10q+r)\cdot\fround{0.1}$ の小数部分を考えて*3、 $$ \begin{aligned} \tfrac9{10} + (q+\tfrac9{10})\cdot 2^{-54} &\ge 1-2^{-4} = \tfrac{15}{16} \\ q+\tfrac9{10} &\ge \tfrac3{80} \cdot 2^{54} \\ &= \tfrac35\cdot 2^{50}. \end{aligned} $$ よって、条件を満たす最小の $q$ は $\ceil{\tfrac35\cdot 2^{50}-\tfrac9{10}} = \tfrac15\cdot(3\cdot2^{50}-2)$ となります。 これにより、最小の $n$ は $10\cdot\tfrac15\cdot(3\cdot 2^{50}-2)+9 = 3\cdot2^{51}+5$ となります。これは、冒頭に挙げた例そのものです。

disclaimer: ここでは、該当の範囲で $(10q+r)\oslash 10 = q$ となることをちゃんと示していません。よしなに。

一般のケース

一般に、$\floor{n\oslash d}$ と $\floor{n\otimes\fround{\tfrac1d}}$ を考えてみます。 $d = (2k+1)\cdot 2^s$ とします。二冪のケースは指数部が変わるだけなので $k\ge 1$ とします。 $$ \tfrac1d = \tfrac{2^t}{2k+1}\times 2^{-(s+t)} $$ と表せます。ここで、$t$ は $2^{p-1}\le\tfrac{2^t}{2k+1}\lt 2^p$ となるように定めます。 すなわち、$t = (p-1)+\ceil{\log_2(2k+1)}$ となります。 $$ \tfrac1d = \underbrace{\tfrac{2^{(p-1)+\ceil{\log_2(2k+1)}}}{2k+1}}_{\mu}\times 2^{-(s+t)}. $$ この $\mu$ の部分を $\floor{\mu}$ または $\ceil{\mu}$ に適切に丸めたものを $m$ として、 $$ \fround{\tfrac1d} = m\times 2^{-(s+t)} $$ となります。 $$ (2^{(p-1)+\ceil{\log_2(2k+1)}})\bmod (2k+1) \le k $$ であれば $m=\floor{\mu}$ に、そうでなければ $m=\ceil{\mu}$ になります。 よって、$0\lt u\le k$ とおき、丸め方向に応じて適切な符号を用いることで $$ \fround{\tfrac1d} = \tfrac{2^t\pm u}{2k+1}\times 2^{-(s+t)} $$ と書くことができます。$2^t$ の持つ素因数が $2$ のみであることと、$2k+1\ge3$ が $2$ 以外の素因数のみを持つことから、$u=0$ にはなりません。

符号ごとに分けて考えます。

まず、$\fround{\tfrac1d} = \tfrac{2^t+u}{2k+1}\times2^{-(s+t)}$ のケースです。$\fround{\tfrac1d}\gt\tfrac1d$ に注意します。 $n = dq + r$ とし、$2^e\le d q+r\lt 2^{e+1}$ とします。

$$ \begin{aligned} (dq+r)\cdot\fround{\tfrac1d} &= (dq+r) \cdot \left(\tfrac{2^t+u}{2k+1}\times 2^{-(s+t)}\right) \\ &= (dq+r) \cdot \tfrac{2^t+u}d\cdot 2^{-t} \\ &= (q+\tfrac rd) \cdot (1+u\cdot 2^{-t}) \\ &= q + \tfrac rd + (q+\tfrac rd)\cdot u\cdot 2^{-t} \\ &= q + \tfrac rd + (q+\tfrac rd)\cdot u\cdot 2^{-(p-1)-\ceil{\log_2(2k+1)}}. \end{aligned} $$ あとは、先ほど同様に小数部分を考えていけばよいでしょう。疲れてきたので割愛です。

次に、$\fround{\tfrac1d} = \tfrac{2^t-u}{2k+1}\times2^{-(s+t)}$ のケースです。$\fround{\tfrac1d}\lt\tfrac1d$ に注意します。

いきなりですが、実はめちゃくちゃ反例があります。下記は Haskell の対話型ツール GHCi での出力です。

ghci> take 10 . map round $ filter (\x -> x * (1.0 / x) /= x / x) [1..]
[49,98,103,107,161,187,196,197,206,214]

$49\otimes\fround{\tfrac1{49}} \lt 1$ ã‚„ $197\otimes\fround{\tfrac1{197}} \lt 1$ などが成り立ちます。先の $\fround{\tfrac1{10}}$ のケースでは身近な範囲（たとえば int に収まる程度）では大丈夫だったりしましたが、こちらに関しては身近な範囲でも全然事故りそうですね。 注意として、浮動小数点数の精度によってこの反例は異なります。下記は Rust の対話型ツール evcxr での出力です*4。

>> (1..).filter(|&x| 1.0 / x as f32 * x as f32 != 1.0).take(10).collect::<Vec<_>>()
[41, 47, 55, 61, 82, 83, 94, 97, 107, 109]
>> (1..).filter(|&x| 1.0 / x as f64 * x as f64 != 1.0).take(10).collect::<Vec<_>>()
[49, 98, 103, 107, 161, 187, 196, 197, 206, 214]

この $49\otimes\fround{\tfrac1{49}}$ のようなものについては、書きたいことがいろいろあるので、別の記事で書くことにします。

また、興味のある読者は floor で満足せず、ceil について考えてみてもよいでしょう。ところで $\ceil{x} = -\floor{-x}$ ã‚„ $\fround{-x} = \fround{x}$ などが成り立ちます。 f32 においては、仮数部が 24 bits しかないため、それなりに十分なテストを行うこともできます。C++ では floating-point promotion によって float が勝手に double になったりするため、Rust などを用いる方が楽そうです。

速度に関して

整数の除算が遅いというのは比較的有名な話かと思います。

そこで、$\sum_{i=1}^{10^{10}} \floor{\tfrac i{10}}$ を、整数型 (i64) と浮動小数点数型 (f64) の $i\oslash 10$ 版と $i\otimes \fround{\tfrac1{10}}$ 版で計測してみました。 また、定数除算最適化の影響を考慮するため、除数を定数として用意するものと実行時に与えるものを比較しました。

もちろん値自体は $O(1)$ で（というか紙とペンで）求めることができて、 $$ \begin{aligned} \sum_{i=1}^{10^{10}} \floor{\tfrac i{10}} &= 9\sum_{i=0}^{10^9-1}i + \sum_{i=1}^{10^9}i \\ &= 10\sum_{i=1}^{10^9-1}i + 10^9 \\ &= 5\cdot (10^9-1)\cdot 10^9 + 10^9 \\ &= 5\cdot 10^{18} - 4\cdot10^9 \\ &= 4999999996000000000 \end{aligned} $$ です。これは $2^{63}$ よりも小さいことを指摘しておきます。

足し合わせる変数には i64 を用い、除算パートのみで上記の違いがある状態にします。 なお、$10^{10} \lt 3\cdot 2^{51}+5$ なので、この範囲で正しい値を計算できることに注意してください。 計測結果は下記の通りでした。

	i64	f64 $\oslash\, {10}$	f64 $\otimes\,\fround{0.1}$
定数	2.26 s	2.86 s	2.24 s
変数	5.75 s	2.90 s	2.20 s

どうやら、いずれも f64 で $\otimes$ を用いるものが最も速かったようです。あらあら。 あと、やはり変数での除算は重いようですね。

今回は floor の方のみ実験しましたが、ceil の方においても適切に丸めた方を用いたりすることで高速に正しい値を得ることが可能そうな気がします。 浮動小数点数型と仲よくなると、適切に証明した上で高速化をすることができそうでうれしいですね。

おきもち

コンテスタント目線では、この手の誤差評価などは “非本質” などと扱われがちな部分な気はしています。 一方で、先ほど言及したように高速化に寄与する場合があり、知っていると役に立つ場面もあるのかなとも思います。

あるいは、テストケースを作成する側に回った際にも、「浮動小数点数型を用いたこうした解法は実は正当である」とか、「誤差がこうなることがわかるので、こういうケースを置いておけば落とせそう」とか、そうした考え方ができるとうれしそうです。AtCoder では（特に double を想定したケースは）考慮されていがちな気はしますよね（long double は見落とされがちな気もしますが）。

ソフトウェア開発の文脈だと、「数学的な証明をしました」で単体テストを書かないのは不安な感じがある気がします。 一方、この手の浮動小数点数が絡むものにおいては、単体テストで適当にいくつかのケースを選んで「大丈夫そうですね」と言うのは不安な気持ちになります。テストしていないケースにおいて誤差でめちゃくちゃになることが容易に起こりうるからです。多少詳しい人でないと、先の例の $\fround{\tfrac1{49}}$ などの反例を挙げることは難しそうですよね。

とはいえ、浮動小数点数に明るくない人（世の中のマジョリティ）にそういう事情を説明してわかってもらうのは大変そうな気がしていて、うまいことやっていく必要がありそうだなという気持ちになったりしています。 この辺の話は正規表現などのトピックでも当てはまりそうですね。

特になにかあったわけではないのですが、競プロ寄りの人々には興味のない話かなと思って薄めの文字で書きました。

あとがき

$\fround{x}$ の記法は、今回思いつきで導入したものになります。 下記の前回の記事では $\hat x$ とか $\operatorname{round}(x)$ とか書いていたものですが、長めの式に対してもよい感じに書けるものがいいなと思ったので考えてみました。

rsk0315.hatenablog.com

経緯としては、最近（文脈自体はあまり関係ないものの）$[\![\texttt{a+b}]\!]$ のような括弧が 使われているの を見て、「括弧をそうやって作るのも “アリ” なんだな」と思ったのがきっかけになります。操作としては round を表すので丸い括弧がいいなと思って、$(\hspace{-.2em}|x|\hspace{-.2em})$ や $\fround{x}$ などを試して決めたものです。 $(\hspace{-.2em}|x|\hspace{-.2em})$ もいいなと思ったのですが、別の文脈*5での使われを思い出したことと、$[\ldots]$ の方が Gauß 記号や $\floor{\ldots}$ $\ceil{\ldots}$ を連想させて、丸めに近いかな？という気持ちがあったことによります。定義は下記です*6。

$\gdef\fround#1{(\hspace{-.2em}[#1]\hspace{-.2em})}$

もしかしたらこの括弧に別の意味があったり、該当の概念に別の記法があったりするのかもしれませんが、あまりサーベイしていません。

おわり

最近は需要のなさそうな記事ばかり書いていますが、何らかの局面で誰かしらの役に立ったらいいなと思っています。あるいは他の人の役に立っていなくても、書いているときの自分が気持ちよくなっているので十分です。

*1:「誤差でぐちゃぐちゃになった意味不明な値」のようなイメージがある人もいるかもしれませんが、数学的に書ける値であることを意識するとよいかもです。

*2:仕様に応じて適切な丸め方向を採用するものとします。

*3:仮数部の偶奇に応じて等号を含むかが変わりますが、一旦等号を含むものから考えてみて、だめならその次のものを採用することにします。

*4:いろいろな対話型ツールがあってたのしいですね。

*5:The SATySFibook の record の部分です。こっちも別の出典があるのかも。分野がやや違いそうなので、同じ記号でも気にしてなくていいと思ってはいます。

*6:論文とかで何らかの独自の記法を使うときは、$\LaTeX$ での書き方を footnote にでも書いておいてくれと思っているタイプです。