为什么对于某些浮点数，a**2 != a * a？

Question

$ python --version
Python 2.7.15

$ type test.py
import random

while True:
    a = random.uniform(0, 1)
    b = a ** 2
    c = a * a
    if b != c:
        print "a = {}".format(a)
        print "a ** 2 = {}".format(b)
        print "a * a = {}".format(c)
        break

$ python test.py
a = 0.145376687586
a ** 2 = 0.0211343812936
a * a = 0.0211343812936

我只能在 Python 的 Windows 构建版本上复现这个问题 - 确切地说：Python 2.7.15 (v2.7.15:ca079a3ea3, Apr 30 2018, 16:30:26) [MSC v.1500 64 bit (AMD64)] on win32。在我安装了 Python 的 Arch Linux 上（Python 2.7.15 (default, May 1 2018, 20:16:04) [GCC 7.3.1 20180406] on linux2），这个循环似乎没有终止，这表明 a**2 = a * a 不变式在那里成立。

这里发生了什么？我知道 IEEE 浮点数有许多误解和特异性（例如，这个并不能回答我的问题），但我不明白规范的哪个部分或者什么样的 ** 实现可以允许这种情况。

针对重复标记的问题：这很可能不是直接的 IEEE 浮点运算问题，更多的是关于 ** 运算符实现的问题。因此，这不是与浮点问题（如精度或结合性）相关的重复问题。

- user1658887

---

使用 repr 显示更多位数，这样我们就可以看到它失败的确切值以及如何失败。 - user2357112

这并不是“浮点数计算是否有问题”的重复。发帖者想知道为什么a ** 2没有被实现为a * a。 - Oliver Charlesworth

@EricPostpischil：乘方不是IEEE 754要求正确舍入的操作。依赖于a**2 == a*a只会带来痛苦。 - user2357112

无论 IEEE 754 是否要求指数运算正确舍入（它确实建议这样做），浮点数的性质 不是 OP 观察到的问题的原因。浮点数操作完全能够支持等式 a**2 == a*a 的恒等性。因此，如果在某些实现中它不能成立，则有其他原因导致，而非浮点数的性质。我并没有建议 OP 或任何人依赖于这个恒等式。我只是陈述了失败的原因不是由于浮点运算的性质。 - Eric Postpischil

欢迎来到浮点数，这里没有什么是精确的。说真的，这种情况经常发生。我们不需要在每个应用程序中都提出一个新的问题。当浮点数以数学建议的方式准确运作时，我们应该更加惊讶。 - jpmc26

Answer 1

Python在浮点数算术方面依赖于底层平台。我假设Python的**运算符使用了一个pow实现（与C中使用的相同）（由user2357112引用Python 2.7.15源代码确认）。

通常，pow部分地使用（近似）对数和指数来实现。这是必要的，因为pow支持非整数参数。（当然，这种一般实现不排除其定义域子集的专业化。）

微软的pow实现是臭名昭著的不好。因此，对于pow(a, 2)，它可能返回一个与a*a不相等的结果。

- Eric Postpischil

---

你可以在这里看到C语言中的pow函数在float.__pow__中的使用：https://github.com/python/cpython/blob/v2.7.15/Objects/floatobject.c#L921。 - user2357112

那个 float_pow 实现当然只适用于 CPython——但是如果你看看 Jython、IronPython 或者 PyPy，你会发现它们在 Java、.NET 和 RPython 中使用了类似的函数。 - abarnert

我有一种预感，在Linux上至少对小整数幂进行了特殊处理，实际上，在a*a，a**2，math.exp(2*math.log(a))中，前两个似乎总是相等的，而第三个则经常偏离。 - Paul Panzer

@PaulPanzer：这并不一定意味着整数幂被特殊处理。在一个完全正确舍入的平台上，exp和log的幂运算应该得到完全相同的结果：math.exp(2 * math.log(a))通常会因为中间值math.log(a)的额外舍入误差而偏离。这个误差可能会被exp调用放大（取决于a的值）。一个好的pow实现应该比exp(y*log(x))做更复杂的事情。 - Mark Dickinson

Answer 2

a ** 2 使用浮点数幂函数（类似于您可以在标准C数学库中找到的函数），能够将任何数字升至任何幂次方。

a * a 只是乘以一次，更适合这种情况，并且不容易出现精度误差（对于整数更为真实），而a ** 2则会出现。

对于浮点数a，如果您想将其升至5次幂，则可以使用以下方法：

a * a * a * a * a

如果重复进行乘法操作，现在容易引起浮点累积错误，并且速度要慢得多，因此你最好使用a**5。

当b较大时，使用a ** b更加高效。但是，由于它使用浮点算法，精度可能会有所不同。

- Jean-François Fabre

---

float.__pow__ 函数并没有使用 math.pow 函数，至少在 CPython 中不是这样的。然而，它们最终都会使用来自 C 标准库的相同 pow 函数（在一堆预检和转换之后，两者并不完全相同），这可能是你应该在这里说的内容。 - abarnert

是的，稍微有些不准确。 - Jean-François Fabre

即使是 * 运算也可能存在精度误差，如果原始数字使用了足够多的位数。2.5 * 2.5 可以得到准确结果，但是 2.123456789 * 2.123456789 将会有一些精度损失。 - Loren Pechtel

@LorenPechtel 2.5*2.5 之所以能够正常运行，是因为它在二进制中可以精确表示，并且只需要很少的位数。尝试运行 2.2*2.2 (4.840000000000001) 就会出现问题。 - jpmc26

@jpmc26 是的。我故意选择了一个浮点数中表示较短的值，以便不会有精度损失。 - Loren Pechtel