0.1 + 0.2が0.3にならないのはなぜか

1. 単精度と倍精度の数は、符号$s$、指数$e$、仮数$f$で1つの数を表す。$p$は先頭の1を含めた仮数の桁数で、単精度は$p=24$、倍精度は$p=53$になる。指数のバイアス$b$は、単精度が127、倍精度が1023である。

   $$x=(-1{)}^s\times \left(1+\frac{f}{2^{p-1}}\right)\times 2^{e-b}$$

2. 0.1は$1.6\times 2^{-4}$なので、単精度の指数は$-4+127=123$。仮数に入りきらない部分を丸めると、単精度の0.1は次の分数に変わる。

   $$\begin{array}{rl}{\mathrm{fl}}_{32}(0.1)& =13421773\times 2^{-27}\\ & =0.100000001490116\dots \end{array}$$

3. 倍精度では53桁で丸めるので、分子が大きくなる。$\mathrm{fl}(x)$は、$x$にいちばん近い倍精度の数を表す。

   $$\mathrm{fl}(0.1)=3602879701896397\times 2^{-55}$$

4. $2^E\le x<2^{E+1}$の範囲では、隣り合う数の間隔はどこでも同じである。$E$は、$x$の2を底とする対数の整数部分にあたる。この間隔をulpと呼ぶ。英語のunit in the last place（最後の桁の1単位）の頭文字だ。

   $$\mathrm{ulp}(x)=2^{\lfloor {\log }_2|x|\rfloor -p+1}$$

5. 間隔を値で割った比は、$2^{-p}$より大きく$2^{1-p}$以下に収まる。いちばん近い数に丸めれば誤差は間隔の半分までなので、相対誤差は$u=2^{-p}$を超えない。

   $$\mathrm{fl}(x)=x(1+\delta ),|\delta |\le 2^{-p}$$

6. 10進で何桁ぶんかは、$p$に${\log }_{10}2$を掛けて求まる。単精度は約7桁、倍精度は約16桁になる。

   $$24{\log }_{10}2=7.22,53{\log }_{10}2=15.95$$

7. 倍精度の0.1と0.2を、丸めずに足す。0.2の分数は0.1の分数のちょうど2倍なので、和は0.1の分数の3倍になる。

   $$\begin{array}{rl}& 3602879701896397\times 3\times 2^{-55}\\ & =10808639105689191\times 2^{-55}\end{array}$$

8. 和は1/4と1/2の間にあるので、倍精度で持てる数の間隔は$2^{-54}$である。$2^{-55}$を単位にした分子は奇数で、和は隣り合う2つの数のちょうど真ん中にある。

   $$\begin{array}{rl}& 10808639105689191\times 2^{-55}\\ & =\left(5404319552844595+\frac{1}{2}\right)\times 2^{-54}\end{array}$$

9. 真ん中に落ちたときは、仮数の最後のビットが0になる側、つまり分子が偶数の側を選ぶ。0.3を倍精度にした数は分子が5404319552844595なので、答えはそれより$2^{-54}$だけ大きい。$\oplus$は計算機の足し算を表す。

   $$\begin{array}{rl}0.1\oplus 0.2& =5404319552844596\times 2^{-54}\\ & =\mathrm{fl}(0.3)+2^{-54}\end{array}$$

単精度で同じ計算をすると、正確な和は$10066329.75\times 2^{-25}$で、真ん中には落ちない。丸めた答えは、0.3を単精度にした数とぴったり重なる。

0.1 + 0.2 == 0.3が成り立つかどうかは、使う精度で変わる。

1. 和を$s$に、丸めで失った分を$c$に持ち、どちらも0から始める。次の項$x_j$から、前の回に失った分を引いておく。

   $$y=x_j-c$$

2. 和に足す。$s$が大きいと、$y$の下の桁は丸めで消える。

   $$t=s+y$$

3. $t-s$は、$y$のうち実際に和に入った部分である。そこから$y$を引くと、消えた部分の符号を変えた値が残る。

   $$c=(t-s)-y$$

4. 和を新しくする。$c$は、次の回の手順1で足し戻される。

   $$s=t$$

5. ゴールドバーグの総説の定理8によれば、こうして求めた和は、各項を$2\epsilon$以内だけ変えた値の正確な和に、$n{\epsilon }^2$の程度の誤差を足したものになる。$\epsilon$は、丸め1回あたりの相対誤差の上限である。

   $$\begin{array}{rl}S& =\sum _{j=1}^n{x}_j(1+{\delta }_j)+O(n{\epsilon }^2)\sum _{j=1}^n|x_j|\\ |{\delta }_j|& \le 2\epsilon \end{array}$$

6. 素直に足したときは、先に足した項ほど多くの丸めを受け、$j$番目の項が$(n-j)\epsilon$近くまで変わりうる。

   $$S=\sum _{j=1}^n{x}_j(1+{\delta }_j),|{\delta }_j|<(n-j)\epsilon$$

単精度のまま$c$を1つ持つだけで、0.1を100万回足した和は100,000になった。1から100万までの和は500,000,489,472で、正確な和500,000,500,000にいちばん近い単精度の数である。

ただし、計算機の足し算をふつうの数の足し算とみなすと、手順3は$c=(s+y-s)-y=0$になってしまう。ゴールドバーグは、そう考えて式を簡単にするコンパイラーの最適化が、この手順を役に立たなくすると注意している。