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

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

原因是您使用的是浮点值，这些值在设计上不是精确或精确 - Python 中的float 刚刚53 个有效位。对于这意味着什么，我建议您阅读文章 What Every Computer Scientist Should Know About Floating-Point Arithmetic或者稍微容易一点的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

应该向下舍入。我猜你的处理器在那里“错误地”舍入了值。

---

对于 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 中，您还可以轻松检查 float 的二进制表示形式:

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

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

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

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

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

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

>>> len('6249999999999997779553950749687116367962969071310727932141320696557418301608777255751192569732666015625')
103