浮点数的计算误差——一个难以忽视的真相（原创）

香酥扇贝

先看一个简单的程序：
Program main
implicit none
real::undef=-9999.9,x=-9999.91
print *,undef,x,undef-x
End program
运行之后结果显示：
-9999.900      -9999.910      9.7656250E-03
与理论值0.01相差了0.002还要多……
可见单精度浮点数只能保证6位有效数字以内误差较小，如果程序计算复杂，涉及到非线性计算，累积误差可能变得很大……

通过debug查看一下undef的数值，发现保存的竟然是-9999.90039，莫名其妙多出来 -0.00039……

简单探讨一下原因：现在的32bit浮点数一般遵从IEEE-754标准，存储方式为最高一位为符号位，然后接着是8位指数，最后23位为尾数。
所以的十进制数都是按照二进制展开的，由于尾数位数有限，因此表达精度有限。
例如：9999.9=2^13+2^10+2^9+2^8+2^3+2^2+2^1+2^0+2^(-1)+2^(-2)+2^(-3)+2^(-6)+2^(-7)+2^(-10)+2^(-11)+...
由于float32存储位数有限，从-11次方开始要截断，因此实际值为
2^13+2^10+2^9+2^8+2^3+2^2+2^1+2^0+2^(-1)+2^(-2)+2^(-3)+2^(-6)+2^(-7)+2^(-10)=9999.8994140625
但还不是9999.90039
再看看临近的稍微大点儿的数：
2^13+2^10+2^9+2^8+2^3+2^2+2^1+2^0+2^(-1)+2^(-2)+2^(-3)+2^(-6)+2^(-7)+2^(-9)=9999.900390625
挺像了吧！比较一下，9999.8994140625与9999.9相差了0.0005859375，
而9999.900390625与9999.9则只相差了0.000390625，更加精确，因此编译器也选择了这个数存储~

如果需要获得更高精度的计算，需要使用64bit双精度，在Linux下编译相同的程序时使用命令：ifort -r8 *.f90
但是双精度比较占用资源，如果只是进行一些简单的计算分析，单精度也够用了~

lqouc

嗯，不错。是楼主自己总结的么？原创的话给加分哈。
话说图挂了。

香酥扇贝

lqouc 发表于 2014-4-16 15:51
嗯，不错。是楼主自己总结的么？原创的话给加分哈。
话说图挂了。

绝对是原创！我总结了大半天了……

香酥扇贝

lqouc 发表于 2014-4-16 15:51
嗯，不错。是楼主自己总结的么？原创的话给加分哈。
话说图挂了。

图去掉了，问题不大，原来的截图就是debugger里面显示的内存中保存的数值~

mofangbao

基本上编程语言都会遇到这个问题 所以直接比较两个值是否相等有时候会出错的~

香酥扇贝

mofangbao 发表于 2014-4-16 15:54
基本上编程语言都会遇到这个问题 所以一般直接比较两个值是否相等有时候会出错的~

是的，要看编译器还有使用的精度是否一致。有时候可以使用|a-b|<eps来代替a==b……

言深深

哈哈，看你的ID很有食欲啊？话说刚刚过去的这个周末吃了不少···

香酥扇贝

言深深 发表于 2014-4-16 16:09
哈哈，看你的ID很有食欲啊？话说刚刚过去的这个周末吃了不少···

哈哈，海鲜很好吃呀！我也爱吃~

言深深

香酥扇贝 发表于 2014-4-16 16:33
哈哈，海鲜很好吃呀！我也爱吃~

哈哈，可不是嘛···可惜大兰州少啊···

渭水听风

过来学习一下