拜读了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位,也就是
eps(0.1)
2^(-56)
ans =
1.3878e-17
ans =
1.3878e-17
y = eps(x) 定义就是计算浮点数集中比
我们知道了间隔,就很容易得到相邻的比x小一点的数,我们记为
因为
所以
同理
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
为了计算
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
因此可以看到
在计算机中
因此
用matlab验证一下
x = 0.1
e = sym(x, 'e')
x =
0.1000
e =
eps/40 + 1/10
与预想结果一致。