对值进行平方的不同方法会导致输出略有不同 - 哪种方法最准确/可靠？

Question

我得到(非常偶尔)C代码和Python代码之间略有不同的计算结果,并设法找到一个例子.在Python中,我得到了这个:

>>> print "%.55f" %\
... (-2.499999999999999555910790149937383830547332763671875 *\
... -2.499999999999999555910790149937383830547332763671875)
6.2499999999999982236431605997495353221893310546875000000
>>> print "%.55f" %\
... ((-2.499999999999999555910790149937383830547332763671875) ** 2)
6.2499999999999973354647408996243029832839965820312500000
>>> print "%.55f" %\
... math.pow(-2.499999999999999555910790149937383830547332763671875, 2)
6.2499999999999973354647408996243029832839965820312500000

而在C中,以下程序:

#include<math.h>
#include<stdio.h>
int main(){
    printf("%.55f\n", -2.499999999999999555910790149937383830547332763671875\
    * -2.499999999999999555910790149937383830547332763671875);
    printf("%.55f\n",\
    pow(-2.499999999999999555910790149937383830547332763671875, 2));
    return 0;
}

给出以下结果:

6.2499999999999982236431605997495353221893310546875000000
6.2499999999999982236431605997495353221893310546875000000

它变得更糟.跑步bc -l,我得到以下内容:

-2.499999999999999555910790149937383830547332763671875 *\
-2.499999999999999555910790149937383830547332763671875 
6.249999999999997779553950749687116367962969071310727

这似乎是正确的结果; 卡西欧的在线高精度计算器表示赞同.

令我担心的,但是,是(var * var),(var ** 2)和math.pow(var, 2)偶尔给的结果稍有不同(我使用Python 2.7.6).

谁知道为什么？

Answer 1

原因是你使用的是浮点值,这在设计上并不准确或精确 - float在Python中只有53个有效位.对于这意味着什么,我建议你阅读文章什么每个计算机科学家应该知道浮点运算或稍微容易每个计算机程序员应该知道浮点.

虽然在Python的情况下,也存在舍入错误:

# incorrect one is:
>>> (6.2499999999999973354647408996243029832839965820312500000).hex()
'0x1.8fffffffffffdp+2'
>>> (6.2499999999999982236431605997495353221893310546875000000).hex()
'0x1.8fffffffffffep+2'

什么时候会有精确的结果

0x1.8fffffffffffd8000000000001p+2

请注意,舍入误差是它离实际舍入点小于2 ^ -100,即:

0x1.8fffffffffffd8000000000000p+2

应该四舍五入,而

0x1.8fffffffffffd7ffffffffffffp+2

应该四舍五入.我猜你的处理器在那里"错误地"舍入了值.

对于C程序,请注意pow(-2.499999999999999555910790149937383830547332763671875, 2),它是一个常量值,它被完全消除,因此程序可以在没有-lm数学库的情况下链接; 我建议你用变量中的值来尝试这个部分:

volatile double a = -2.499999999999999555910790149937383830547332763671875;
printf("%.55f", pow(a, 2));

这可能会提供不同的结果.

但Python也有一个任意精度的十进制库,decimal.Decimal:

with localcontext() as ctx:
    # set the precision to 200 significant digits
    ctx.prec = 200
    result = decimal.Decimal(
       '2.499999999999999555910790149937383830547332763671875') ** 2)

print(result)

版画

6.24999999999999777955395074968711636796296907131072793214132069
6557418301608777255751192569732666015625

这是正确的结果,而不是你得到的结果bc; 在我的计算机bc -l上将比例设置为20位有效数字,但它可以通过scale变量进行调整:

% bc
bc 1.06.95
Copyright 1991-1994, 1997, 1998, 2000, 2004, 2006 Free Software Foundation, Inc.
This is free software with ABSOLUTELY NO WARRANTY.
For details type `warranty'. 
scale = 200 
2.499999999999999555910790149937383830547332763671875 * 2.499999999999999555910790149937383830547332763671875
6.249999999999997779553950749687116367962969071310727932141320696557\
418301608777255751192569732666015625

结果与Python相同.(注意100的精度不够,因为结果值有104个有效十进制数字).

您也可以通过python检查数字,没有任何导入 - Python也有任意精度整数运算:

>>> 2499999999999999555910790149937383830547332763671875 * \
... 2499999999999999555910790149937383830547332763671875
624999999999999777955395074968711636796296907131072793214132069
6557418301608777255751192569732666015625L

在Python 3中,您还可以轻松检查浮点数的二进制表示:

>>> (2.499999999999999555910790149937383830547332763671875).hex()
'0x1.3ffffffffffffp+1'

此数字的平方应与具有相同数字平方的整数具有相同的数字:

>>> hex(0x13ffffffffffff * 0x13ffffffffffff)
'0x18fffffffffffd8000000000001'

但很明显,这仅仅在53个有效位中无法表示; 相反,它需要两倍的位才能完全表示:

>>> log2(0x13ffffffffffff * 0x13ffffffffffff)
104.64385618977472

顺便提一下,由于fs,这几乎是十进制表示所需的非零数字的数量:

>>> len('6249999999999997779553950749687116367962969071310727932141320696557418301608777255751192569732666015625')
103