0.1 + 0.9がちょうど1.0になるやつ、試してみた

# こんなツイートが流れてきました

<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">今すごいびっくりしてしまったんですが、<br> double a = 0.1, b = 0.9;<br>に対して、<br> assert(a + b == 1.0);<br>がC++で成立するの非自明じゃないですか？なんで？？？</p>&mdash; はむこ (@hamko_intel) <a href="https://twitter.com/hamko_intel/status/1231933387209756674?ref_src=twsrc%5Etfw">February 24, 2020</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script> 

嘘やん。ならないって思ってた。

## そもそも何の話か分からない人へ
「え? 0.1 + 0.9がちょうど1.0になるのは当たり前では?」って方もいるかもしれません。
ざっくり言うと、0.1も0.9も、2進数の小数で表すと無限に続く循環小数になって、double型などで表すと途中で桁を打ち切ることになるので、足してちょうど1.0になるのは不思議だなぁ、という話をしています。

この記事では、浮動小数点数に誤差が出る理由は解説しませんので、知りたい方は「浮動小数点数 誤差」などで検索してください。例えば[この記事](https://qiita.com/angel_p_57/items/24078ba4aa5881805ab2)がおすすめです。

# やってみよう
まず、`double`と`uint64_t`のunionを作って、`double`のビット列を取り出して、符号部、指数部、仮数部を取り出します。[^1]
次に、小さい方(0.1)の仮数をシフトして、大きい方(0.9)に合わせてから、仮数を加算しました。[^2]


[^1]: Wikipediaの[IEEE 754](https://ja.wikipedia.org/wiki/IEEE_754)を参考にしました。IEEE 754のビット列をどういうエンディアンで詰めるとかは、情報を見つけられなかったため、私の環境(x64)に合うよう実装しています。

[^2]: 浮動小数点数の加算は[ディジタル回路設計とコンピュータアーキテクチャ 第2版](https://www.shoeisha.co.jp/book/detail/9784798147529)の電子書籍の5.3節にある方法を参考にしています。

```cpp:みんな大好きC++
#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になる様子が見て取れます。


# ところで

<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">自宅のWindowsのCLion(MinGW)でやってみたらAssertion failedになってしまいました。<br> double a = 0.1, b = 0.9;<br> cout &lt;&lt; (a + b - 1.0) &lt;&lt; endl;<br>こうすると 2.77556e-017 が出力されました。<br><br>Wandboxでは等式が成立したので、この浮動小数点演算は処理系に依存している気がします。</p>&mdash; Lillian (@Lily0727K) <a href="https://twitter.com/Lily0727K/status/1231944416216014848?ref_src=twsrc%5Etfw">February 24, 2020</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script> 

<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">確かに処理系依存には違いないです。x86_64 Linux で -mfpmath=387 をつけて g++ でコンパイルすれば同様の結果になったので、x87 の 80bit floating point によるものだと思われます。</p>&mdash; ぐらふぃ (@grafi_tt) <a href="https://twitter.com/grafi_tt/status/1231957093021700096?ref_src=twsrc%5Etfw">February 24, 2020</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script> 

確かに、g++に`-mfpmath=387`を付けてコンパイルすると、assertが通らなくなりました。

```cpp:問題のコード
#include <cassert>

int main()
{
    double a = 0.1, b = 0.9;
    assert(a + b == 1.0);
}
```

また、次のように`a + b`を一旦変数に入れると、`-mfpmath=387`を付けても、assertは通るようになりました。

```cpp:変数に受けてみる
#include <cassert>

int main()
{
    double a = 0.1, b = 0.9;
    double c = a + b;
    assert(c == 1.0);
}
```

80 bit floating pointについて、このあたりを参考に考察しました。

- [浮動小数点数演算に関するメモ](http://syohex.hatenablog.com/entry/20110831/1314802064)
- [浮動小数点演算ではまった話](http://0xcc.net/blog/archives/000164.html)
- [x86における浮動小数点演算の精度の制御と、80bit 浮動小数点演算問題](http://www.takeoka.org/~take/cpu/num/float80bit.html)

そうすると、さっきは、小さい方(a = 0.1)の仮数を右シフトして、桁落ちさせてから考えましたが、内部には十分な精度の演算器があるので、小さい方の仮数を桁落ちさせる必要がないと考えられます。
なので、小さい方の仮数を右シフトするかわりに、大きい方(b = 0.9)の仮数を左シフトして計算してみました。

```cpp:main関数にこれ付け足す
    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浮動小数点演算ではそうならない理由を見ていきました。