浮点数精确运算的分析和解决办法

浮点数精确运算的分析和解决办法

    1.01 + 2.01 = 3.02
    2.01 * 2.01 = 4 0401
   不知你注意没有，这个很寻常的等式，你如果将它放在C++中，Java中，Basic中，它
居然是不成立的。计算机在开玩笑吗？噢，对了，隐约记得这好象是浮点数的问题，似乎
很多很多年前，老师说过。还有某位姓林的先生在某本书里提过=0的判断。
   嗯，如果你不遇到此问题，那你完全可以把它抛到火星上去，可惜，偶不好彩，这样的
问题，被俺遇到了。唉！
   why？how？
  
   没办法，硬着头皮，从头开始。

一：为何不成立？Why？
   这得从浮点数的在计算机内的存储开始说起，我这里闲话少说。我们只谈双精度double
数(至于float，基本上是五十步和一百步的区别)。
   双精度数在计算机内的表示方式是：(三部分组成)
   符号(正或负)  阶码(2的N次幂)   尾数(大于等于1小于2的数)
   比如： -(符号) 1.01(尾数) * 2~1(N = 1)  = - 2.02
   具体到计算机的存储单元：双精度数共占8字节(64bit)
   符号位(占1个bit) 阶码(11个bit) 尾数(52个bit)
解释一下：
   符号位：0表示正 1表示负
   阶码：是一个偏移量，1023的偏移量,它的1023相当于0,小于1023时为负,
大于1023时为正,如:10000000001表示指数为1025 - 1023 = 2,表示真值为2^2。

好了，知道了原理，我们开始分析上述等式为何为不等。
(相应数的存储值，可以简单用C语言的指针方式取出)
1.01 表示为：
0   0111111 1111    0000 00101000 11110101 11000010 10001111 01011100 00101001 
2.01 表示为；
0   1000000 0000    0000 00010100 01111010 11100001 01000111 10101110 00010100
3.02 表示为：
0   1000000 0000    1000 00101000 11110101 11000010 10001111 01011100 00101001
2.01+1.01 在编程语言中的计算结果 表示为：
0   1000000 0000    1000 00101000 11110101 11000010 10001111 01011100 00101000
好了，我们可以比较一下3.02和计算结果，果然有所不同，只不过最后一个bit不同嘿。
为了验证一下，可以用手工计算一下2.01+1.01：
先把1.01的幂次变为1(与2.01的阶码相同)，于是，将尾数右移一位。得到：
  1000 00010100 01111010 11100001 01000111 10101110 000101001
加上2.01的尾数。
  0000 00010100 01111010 11100001 01000111 10101110 00010100
得到：
  1000 00101000 11110101 11000010 10001111 01011100 00101000
嗯，与计算机的计算结果相同，我们的运算思路是正确的。

因此，结论出来了，因为浮点数在计算机内的存储存在偏差，导致运算时，与实际期望的结
果不同。很多时候，你可以不理它，但是，可以肯定负责任的说，发射卫星的运算时，你
需要知道，否则，卫星一转眼就不见了。

二：不成立的的原因找到了，那怎么解决这个问题呢，How？
一个简单的解决办法是：
不要用浮点数来存储浮点，对于VC,Java,Basic，最好的办法是用Decimal来保存它。
下面是分别的实现：(以加法为例，其它四则运算处理相同)
VC中：
double doublAdd(double dbl1, double dbl2)
{
 double dblResult;
 DECIMAL dec1,dec2,decResult;
 ::VarDecFromR8(dbl1,&dec1);
 ::VarDecFromR8(dbl2,&dec2);
 ::VarDecAdd(&dec1,&dec2,&decResult);
 ::VarR8FromDec(&decResult,&dblResult);
 return dblResult;
}
VB中：
Private Function doubleAdd(ByVal dbl1 As Double, ByVal dbl2 As Double) As Double
    doubleAdd = CDec(dbl1) + CDec(dbl2)
End Function

Java中：
public static double add(double v1, double v2) {
    BigDecimal b1 = new BigDecimal(Double.toString(v1));
    BigDecimal b2 = new BigDecimal(Double.toString(v2));
    return b1.add(b2).doubleValue();
}
解决思路就是：用其它精确的表示法来存储浮点数，就这么简单。
注意：VC示例中，VarDecFromR8是做了手脚地，如果能直接用VarDecFromStr那更好。

三：在C/C++中，似乎很不情愿看到类似上例中的代码，因为它看起来很低效，还有其它方法吗？
好象还有，对了，只是好象。
我们再来看看双精度数的表示法：
尾数一共有52个bit，也就是最小能表示的数是 2^-52，取对数可得出，约是
在小数点后16位，那也就是说小数点后15位是可以精确表示的，加上前置的默认1，一共有16位
数字是精确可靠的。
我们来试验一下，看上述结论是否成立。
看看VC调试器的显示值。
2.01 的显示值： 2.0099999999999998
如果只取16位有效数字，那么将最后一位8四舍五入，我们得到正确的表示。
好了，这能说明什么呢？

四：我们先看比较简单的加，减法运算。
对于加法：dbl1 + dbl2：
假设dbl1=1.01 那么，16减去整数位1，我们可以假定，在计算机表示中：
小数点后的15位都是精确的。
假设dbl2=100.01 那么 16-3，假定小数点后13位是精确的。
凭经验我们可以知道，两个小数相加，小数点后的精度不会大于精度销大的一个。
所以，我们判定得出结果的精确度可以用较大的一个为准。
于是，将得出的结果，去掉不精确的位数，则应该可以得到准确值。
VC实现如下：
#define DELTA_RATE  16
int getRound(double dbl)
{
 COleVariant var(dbl);
 COleVariant varForLog(dbl);
 ::VarRound(&varForLog,0,&varForLog);
 int nIntCount = log10(varForLog.dblVal>0?varForLog.dblVal:-varForLog.dblVal) + 1;
 int nRound = DELTA_RATE - nIntCount;
 return nRound;
}
double doublAdd2(double dbl1, double dbl2)
{
 COleVariant var(dbl1+dbl2);
 int r1 = getRound(dbl1);
 int r2 = getRound(dbl2);
 ::VarRound(&var,max(r1,r2),&var);
 return var.dblVal;
}
做过一些实验，好象是正确的。同理可以实现doubleSub2的函数。
注意：这里并不用下面五所提的取精度的方式，因为取精度的运算更低效。

五：对于乘除法呢？问题有些复杂，先找出一个需要处理的例子。
如：2.01*2.01=4.0401。
试了一下，不成立。
用方法一的Decimal方式测试，可以通过。
那么方法二呢？
再做假设吧，假设dbl1有两位小数，dbl2也有两位小数，按理论，
可得出相乘后，最大可能是2+2位小数。那么，我们按照 4位小数
进行Round处理，可能会得出正确的结果。
实际上，要取一个双精度的10进制表达的小数位，我没有找到什么好办法，
我能想到的：也就是将数字转为字串，然后查找.后的位数。这样，显然是
非常低效的，这里，我就不再写出代码了。

六：比较方法一和方法二。方法二并不高效，并且还有一些不定因素，所以，
最好采用方法一来统一处理浮点数的运算。
至于效率，实际上最佳方法是从程序的设计着手，将double从程序中去除掉。
比如在VC中，可以用Variant::Decimal来彻底替换double，这样，就不存在
中间的转换了，效率自然就提高了。有关Decimal的常用函数是：
VarDecFromStr VarDecAdd VarDecSub VarDecMul VarDecDiv ……
VarBstrFromDec
至于Java和VB，也可以方便的找到相应函数。

很想找到一种更好的方法，总觉得用Decimal来进行运算很不爽，但真的没找到？
其实呢，做了一下测试，Decimal的运算并不慢，如果可以将内部存储改为Decimal，
那就可以彻底解决问题了。