こんなツイートが流れてきました
今すごいびっくりしてしまったんですが、
— はむこ (@hamko_intel) February 24, 2020
double a = 0.1, b = 0.9;
に対して、
assert(a + b == 1.0);
がC++で成立するの非自明じゃないですか？なんで？？？
嘘やん。ならないって思ってた。
そもそも何の話か分からない人へ
「え? 0.1 + 0.9がちょうど1.0になるのは当たり前では?」って方もいるかもしれません。
ざっくり言うと、0.1も0.9も、2進数の小数で表すと無限に続く循環小数になって、double型などで表すと途中で桁を打ち切ることになるので、足してちょうど1.0になるのは不思議だなぁ、という話をしています。
この記事では、浮動小数点数に誤差が出る理由は解説しませんので、知りたい方は「浮動小数点数 誤差」などで検索してください。例えばこの記事がおすすめです。
やってみよう
まず、doubleとuint64_tのunionを作って、doubleのビット列を取り出して、符号部、指数部、仮数部を取り出します。1
次に、小さい方(0.1)の仮数をシフトして、大きい方(0.9)に合わせてから、仮数を加算しました。2
# include <iostream>
# include <string>
# include <cassert>
std::string bits(uint64_t n, int bit_length)
{
std::string s(bit_length, '0');
for (auto i = bit_length - 1; i >= 0; i--) {
if (n & 1) {
s[i] = '1';
}
n >>= 1;
}
return s;
}
union Double {
double d;
uint64_t u;
Double(double d) : d(d) {}
int sign() {
return (u >> 63) ? -1 : 1;
}
char sign_char() {
return sign() == 1 ? '+' : '-';
}
int exponent() {
return (static_cast<int>(u >> 52) & ((1 << 11) - 1)) - ((1 << 10) - 1);
}
uint64_t fraction() {
return u & ((1ull << 52) - 1ull);
}
uint64_t fraction_noneconomized() {
return fraction() | (1ull << 52);
}
std::string fraction_bits() {
return bits(fraction(), 52);
}
};
int main()
{
using std::cout;
using std::endl;
Double d0(0.1), d1(0.9);
cout << d0.d << '\t' << d0.sign_char() << "1." << d0.fraction_bits() << " x 2^" << d0.exponent() << endl;
cout << d1.d << '\t' << d1.sign_char() << "1." << d1.fraction_bits() << " x 2^" << d1.exponent() << endl;
assert(d0.d > 0. && d0.d < d1.d);
auto exponent_diff = d1.exponent() - d0.exponent();
auto d0_bits = d0.fraction_noneconomized() >> exponent_diff;
auto d1_bits = d1.fraction_noneconomized();
cout << endl;
cout << " " << bits(d0_bits, 53) << " x 2^" << (d1.exponent() - 52) << endl;
cout << "+ " << bits(d1_bits, 53) << " x 2^" << (d1.exponent() - 52) << endl;
cout << std::string(63, '-') << endl;
cout << " " << bits(d0_bits + d1_bits, 54) << " x 2^" << (d1.exponent() - 52) << endl;
}
実行結果
0.1 +1.1001100110011001100110011001100110011001100110011010 x 2^-4
0.9 +1.1100110011001100110011001100110011001100110011001101 x 2^-1
00011001100110011001100110011001100110011001100110011 x 2^-53
+ 11100110011001100110011001100110011001100110011001101 x 2^-53
---------------------------------------------------------------
100000000000000000000000000000000000000000000000000000 x 2^-53
0.1の仮数は1001 1001 ...の循環小数になっていて、最後が1010になっています。これは、本来は
... 1001 1001 | 1001 ... (縦線は仮数部の桁数の切れ目を表す)
と続きますが、桁数が切れるところで四捨五入(2進数なので0捨1入)された結果、切り上げられたと考えられます。
0.9の仮数は1100 1100 ...の循環小数になっていますが、これも、最後に0捨1入で切り上げられたと考えられます。
足し算は、指数を大きい方(0.1)に合わせるよう仮数部をシフトし、仮数部の足し算をします。
とてもとても綺麗に、ちょうど1になる様子が見て取れます。
ところで
自宅のWindowsのCLion(MinGW)でやってみたらAssertion failedになってしまいました。
— Lillian (@Lily0727K) February 24, 2020
double a = 0.1, b = 0.9;
cout << (a + b - 1.0) << endl;
こうすると 2.77556e-017 が出力されました。
Wandboxでは等式が成立したので、この浮動小数点演算は処理系に依存している気がします。
確かに処理系依存には違いないです。x86_64 Linux で -mfpmath=387 をつけて g++ でコンパイルすれば同様の結果になったので、x87 の 80bit floating point によるものだと思われます。
— ぐらふぃ (@grafi_tt) February 24, 2020
確かに、g++に-mfpmath=387を付けてコンパイルすると、assertが通らなくなりました。
# include <cassert>
int main()
{
double a = 0.1, b = 0.9;
assert(a + b == 1.0);
}
また、次のようにa + bを一旦変数に入れると、-mfpmath=387を付けても、assertは通るようになりました。
# include <cassert>
int main()
{
double a = 0.1, b = 0.9;
double c = a + b;
assert(c == 1.0);
}
80 bit floating pointについて、このあたりを参考に考察しました。
そうすると、さっきは、小さい方(a = 0.1)の仮数を右シフトして、桁落ちさせてから考えましたが、内部には十分な精度の演算器があるので、小さい方の仮数を桁落ちさせる必要がないと考えられます。
なので、小さい方の仮数を右シフトするかわりに、大きい方(b = 0.9)の仮数を左シフトして計算してみました。
auto d0_bits_ext = d0.fraction_noneconomized();
auto d1_bits_ext = d1.fraction_noneconomized() << exponent_diff;
cout << endl;
cout << " " << bits(d0_bits_ext, 63) << endl;
cout << "+ " << bits(d1_bits_ext, 63) << endl;
cout << std::string(65, '-') << endl;
cout << " " << bits(d0_bits_ext + d1_bits_ext, 64) << endl;
これを実行すると
000000000011001100110011001100110011001100110011001100110011010
+ 000000011100110011001100110011001100110011001100110011001101000
-----------------------------------------------------------------
0000000100000000000000000000000000000000000000000000000000000010
と、右側に不穏なものが残ります。
なので、ちょうど1.0にはならずassert(a + b == 1.0);は通らなかったと考えられます。
一方で、ビットを数えると、不穏な1はdoubleの仮数の中には収まりません。
なので、c = a + b; assert(c == 1.0);は通ったと考えられます。
まとめ
0.1 + 0.9がちょうど1.0になる理由、また、80bit浮動小数点演算ではそうならない理由を見ていきました。
-
WikipediaのIEEE 754を参考にしました。IEEE 754のビット列をどういうエンディアンで詰めるとかは、情報を見つけられなかったため、私の環境(x64)に合うよう実装しています。 ↩
-
浮動小数点数の加算はディジタル回路設計とコンピュータアーキテクチャ 第2版の電子書籍の5.3節にある方法を参考にしています。 ↩