关于浮点数的误差理解

拜读了Cleve Moler写的博客《Bank Format and Metric Socket Wrenches》，对浮点数又有了进一步的理解。原文的链接地址为
http://blogs.mathworks.com/cleve/2017/04/10/bank-format-and-metric-socket-wrenches/
本文的内容是对这篇文章部分翻译以及自己的一些理解，让大家能够更加容易看懂。

从format bank这个问题而来

我们知道matlab中有一些命令可以控制显示格式，比如：
format short： 默认，显示较少的有效数字
format long：显示较多的有效数字
format hex：用16进制显示
format short g：用8进制显示
这些都比较常用，在金融银行业中用的比较多的就是format bank，这个命令保留小数点后2位（四舍五入，保留2位小数）

注意上面的命令只是用来控制显示的格式，而不影响计算的精度

在使用format bank时可能会遇到让你觉得奇怪的现象，如果你对计算机内的浮点数表示不了解的话。以下是一个奇怪的例子

format long
x = (5:10:45)' / 1000
y = 1 + x
z = 2 + x

x =

0.005000000000000
0.015000000000000
0.025000000000000
0.035000000000000
0.045000000000000

y =

1.005000000000000
1.015000000000000
1.025000000000000
1.035000000000000
1.045000000000000

z =

2.005000000000000
2.015000000000000
2.025000000000000
2.035000000000000
2.045000000000000

format bank
x
y 
z

x =
0.01
0.01
0.03
0.04
0.04

y =
1.00
1.01
1.02
1.03
1.04

z =
2.00
2.02
2.02
2.04
2.04

是不是觉很奇怪，变量x, y, z中有些数字四舍五入了，有些没有四舍五入。其实这就是10进制小数不能完全用2进制浮点数精确表示引起的。

虽然在matlab中显示的数字等于 0.015000000000，好像很精确了，但其实在matlab内部这个数却还没有达到0.015，比0.015略小一点。 下面可以证实这一点

exact_x = sym(x, 'e')
exact_y = sym(y, 'e')
exact_z = sym(z, 'e')

sym(x, ‘e’) 返回的就是在matlab中表示的数值x，结果的类型是一个符号数。

exact_x =
eps/2133 + 1/200
3/200 - eps/400
eps/160 + 1/40
(3*eps)/200 + 7/200
9/200 - (3*eps)/400

exact_y =

201/200 - (12*eps)/25
203/200 - (11*eps)/25
41/40 - (2*eps)/5
207/200 - (9*eps)/25
209/200 - (8*eps)/25

exact_z =

401/200 - (12*eps)/25
(14*eps)/25 + 403/200
81/40 - (2*eps)/5
(16*eps)/25 + 407/200
409/200 - (8*eps)/25

可以看到所有的数都比理论值多一点或者少一点，其实这个误差是怎么来的，为什么会有这个误差，其实就跟浮点数的表示有关。我很早写过一篇关于浮点数如何表示的文章，大家可以参看一下，文章地址http://blog.csdn.net/whoispo/article/details/6312782。在下一节中，我会计算这个误差，大家最好能够理解浮点数的表示。

因此对于上面的变量x, y, z，有些进行了四舍五入，有些没有的原因，就比较好理解了。数字比0.5稍小的，就舍去了，比0.5稍大的就入位。

format short 
signx = sign(double(exact_x - x))
signy = sign(double(exact_y - y))
signz = sign(double(exact_z - z))

signx =
1
-1
1
1
-1

signy =
-1
-1
-1
-1
-1

signz =
-1
1
-1
1
-1

等于1就表示matlab中表示值比理论值大，因此入位，等于-1就表示matlab的表示值比理论值小，就舍去。大家看看是不是和format bank显示的结果一致。

计算浮点数的误差

下面就来计算一下浮点数 0.1 的误差。

我们知道0.1在1/16和1/8之间，因此可以确定0.1的浮点数的指数部分应该就是1/16，剩下的就是尾数部分。一般而言，如果不对matlab的类型做特殊说明，matlab的类型大多都是double类型。double类型的尾数部分有52位，也就是$2^{-52}$，再加上指数部分的$2^{-4}$，因此对于1/16到1/8之间的数，相邻两个数之间的间隔为$2^{-56}$，也即0.1附近的数的间隔为$2^{-56}$。用matlab验证一下

eps(0.1)
2^(-56)

ans =

1.3878e-17

ans =

1.3878e-17

y = eps(x) 定义就是计算浮点数集中比$x^*$大一点的数与$x^*$的差值，也就是计算x附近的浮点数的间隔。注意我在这里的x右上角加了一个”*”，表示计算机中x的值。

我们知道了间隔，就很容易得到相邻的比x小一点的数，我们记为$x^-$，和相零的比x大一点的数了，我们记为$x^+$。$x^*$取哪一个数，取决于x离$x^-$近还是离$x^+$近。

$x^-$和$x^+$的计算

$$x^- = floor(\frac{0.1-1/16}{gap})*gap+1/16$$
因为
$$floor(\frac{-1/16}{gap}) * gap + 1/16 = 0$$
所以
$$x^- = floor(\frac{0.1}{gap})*gap$$
同理
$$x^+ = ceil(\frac{0.1}{gap})*gap$$

x = 2^(-56)
symx = sym(x)
x0 = floor(0.1/symx)*symx
x1 = ceil(0.1/symx) * symx

x =

1.3878e-17

symx =

1/72057594037927936

x0 =

7205759403792793/72057594037927936

x1 =

3602879701896397/36028797018963968

为了计算$x$离$x^-$近还是离$x^+$近，下面只能采用matlab的符号计算，为了避免丢失精度，因为x0和x1已经是符号变量了，所以可以直接构造符号变量c。

c = [x0 0.1 x1]'

c =

7205759403792793/72057594037927936
1/10
3602879701896397/36028797018963968

这里的c(2)精确等于0.1，因为这是一个符号。用vpa函数计算精确值，这是对符号数进行高精度计算的方法，可以任意指定精度

vpa(c, 50)

ans =

0.099999999999999991673327315311325946822762489318848
0.1
0.1000000000000000055511151231257827021181583404541

从这里可以看到，用双精度浮点数确实不能精确表示0.1

double(c(1))
double(c(3))

ans =
0.1000

ans =
0.1000

x = 2^(-56);
x = sym(x);
e1 = (0.1 - c(1))/x;
e2 = (c(3) - 0.1)/x;
vpa(e1)
vpa(e2)

ans =

0.6

ans =

0.4

因此可以看到$x$还是离$x^+$近一点。所以0.1在计算机中的存储的值为

$$x^* = 1/10 + eps(0.1) * 0.4$$

在计算机中

$$eps = eps(1.0) = 2^{-52}$$
$$eps(0.1) = 2^{-56}$$

因此

$$x^* = 1/10 + eps / 40$$

用matlab验证一下

x = 0.1
e = sym(x, 'e')

x =
0.1000

e =

eps/40 + 1/10

与预想结果一致。