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

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

这似乎是正确的结果； Casio's online high precision calculator 同意。

然而，让我担心的是 (var * var)、(var ** 2) 和 math.pow(var, 2) 偶尔会给出略微不同的结果（我使用的是 Python 2.7.6）。

有人知道为什么吗？

Answer 1

原因是您使用的是浮点值，这些值在设计上不是精确或精确 - float 在 Python 中只有 53 个有效位。对于这意味着什么，我建议您阅读文章What Every Computer Scientist Should Know About Floating-Point Arithmetic or the slightly easier What every computer programmer should know about floating point。

虽然在 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

应向下舍入。我猜你的 processor 舍入值 "incorrectly" 那里。

对于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

在Python3中你也可以轻松查看浮点数的二进制表示：

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

这个数字的平方应该与具有相同数字平方的整数具有相同的数字：

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

但显然这不能仅用 53 个有效位来精确表示；相反，它需要两倍的位才能准确表示：

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

顺便说一句，由于 fs，这几乎与十进制表示所需的非零数字数量相同：

>>> len('6249999999999997779553950749687116367962969071310727932141320696557418301608777255751192569732666015625')
103

Answer 2

@AnttiHaapala 解释了为什么需要进行舍入 - 结果需要超过 53 位才能获得精确结果。理想情况下，所有三个版本的计算都会得到相同的结果。它们应该全部舍入到相同的值。但并非所有操作系统和编译器都能计算出正确的舍入结果。

我碰巧得到了所有三个版本的相同答案。

Python 2.7.8 (default, Oct 20 2014, 15:05:19) 
[GCC 4.9.1] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import math
>>> a=-2.499999999999999555910790149937383830547332763671875
>>> print "%.55f" % (a*a)
6.2499999999999982236431605997495353221893310546875000000
>>> print "%.55f" % (a**2)
6.2499999999999982236431605997495353221893310546875000000
>>> print "%.55f" % (math.pow(a,2))
6.2499999999999982236431605997495353221893310546875000000

Python 只是使用操作系统和编译器提供的。对于二进制浮点计算，MPFR 库确实保证正确舍入结果。 GCC 使用 MPFR 在编译时评估常量表达式。

我维护 gmpy2 扩展模块，它提供从 Python 访问 MPFR 库的权限。

Answer 3

谢谢你们的回答，伙计们，感谢他们，我更深入地挖掘，并设法更好地理解发生了什么。 @casevh 是正确的，问题不在于 Python，而在于底层 C 库。

我正在使用 Mac，并安装了 gcc 5.0，我用它来编译我提供的示例 C 代码：

$ gcc --version
gcc (GCC) 5.0.0 20141005 (experimental)
Copyright (C) 2014 Free Software Foundation, Inc.

但是 python 解释器正在使用不同的 C 编译器版本：

Python 2.7.8 (v2.7.8:ee879c0ffa11, Jun 29 2014, 21:07:35) 
[GCC 4.2.1 (Apple Inc. build 5666) (dot 3)] on darwin

Apple 提供的 C 编译器是另一个版本：

$ cc --version
Apple LLVM version 6.0 (clang-600.0.56) (based on LLVM 3.5svn)
Target: x86_64-apple-darwin14.1.0
Thread model: posix

出于兴趣，我使用 cc 编译了相同的 C 代码，输出为：

6.2499999999999982236431605997495353221893310546875000000
6.2499999999999973354647408996243029832839965820312500000

所以问题出在 Macs 提供的 C 库版本上。 Python 似乎将 C 数学库中的 pow 函数用于 v ** 2 和 math.pow(v, 2)。

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

Different methods to square a value result in slightly different output - which method is most accurate / reliable?

python

math

python-2.7