matlab - (0.3)^3 == (0.3)*(0.3)*(0.3) 在matlab中返回false？

Question

我试图了解 MATLAB 中基本算术运算的舍入误差，我遇到了以下奇怪的例子。

(0.3)^3 == (0.3)*(0.3)*(0.3)

ans = 0

我想确切地知道左侧是如何计算的。MATLAB 文档建议对于整数幂，使用“平方求幂”算法。

“矩阵幂。X^p 是 X 的幂 p，如果 p 是标量。如果 p 是整数，则幂通过重复平方计算。”

所以我假设(0.3)^3并且(0.3)*(0.3)^2会返回相同的值。但这种情况并非如此。如何解释舍入误差的差异？

Answer 1

我对 MATLAB 一无所知，但我在 Ruby 中尝试过：

irb> 0.3 ** 3
  => 0.026999999999999996
irb> 0.3 * 0.3 * 0.3
  => 0.027

根据Ruby 源代码，如果左侧操作数是浮点数，则求幂运算符将右侧操作数转换为浮点数，然后调用标准 C 函数pow()。该函数的float变体pow()必须实现更复杂的算法来处理非整数指数，这将使用导致舍入误差的操作。也许 MATLAB 的工作原理类似。

Answer 2

有趣的是，标量^似乎是使用实现的，pow而矩阵^是使用平方和乘法实现的。以机智：

octave:13> format hex
octave:14> 0.3^3
ans = 3f9ba5e353f7ced8
octave:15> 0.3*0.3*0.3
ans = 3f9ba5e353f7ced9
octave:20> [0.3 0;0 0.3]^3
ans =

  3f9ba5e353f7ced9  0000000000000000
  0000000000000000  3f9ba5e353f7ced9

octave:21> [0.3 0;0 0.3] * [0.3 0;0 0.3] * [0.3 0;0 0.3]
ans =

  3f9ba5e353f7ced9  0000000000000000
  0000000000000000  3f9ba5e353f7ced9

这可以通过在 gdb 下运行 octave 并在pow.

在 matlab 中可能也是如此，但我无法真正验证。

Answer 3

这是一个小测试程序，它遵循pow()来自Source/Intel/xmm_power.cApple 的系统Libm-2026在这种情况下所做的事情：

#include <stdio.h>
int main() {
    // basically lines 1130-1157 of xmm_power.c, modified a bit to remove
    // irrelevant things

    double x = .3;
    int i = 3;

    //calculate ix = f**i
    long double ix = 1.0, lx = (long double) x;

    //calculate x**i by doing lots of multiplication
    int mask = 1;

    //for each of the bits set in i, multiply ix by x**(2**bit_position)
    while(i != 0)
    {
        if( i & mask )
        {
            ix *= lx;
            i -= mask;
        }
        mask += mask;
        lx *= lx; // In double this might overflow spuriously, but not in long double
    }

    printf("%.40f\n", (double) ix);
}

这打印出来，这与我在 Matlab 和Python 中0.0269999999999999962252417162744677625597得到的结果一致（我们知道后者只是调用这个代码）。相比之下，对我来说 gets ，如果你只是要求打印到那么多小数位，这与你得到的结果是一样的，所以大概是最接近的双精度数。.3 ^ 3.3 ** 3.3 * .3 * .30.02699999999999999969468866822808195138350.027

于是就有了算法。我们可以准确地跟踪每一步设置的值，但如果使用不同的算法，它会舍入到一个非常小的数字也就不足为奇了。

Answer 4

感谢@Dougal，我发现了这个：

#include <stdio.h>
int main() {
  double x = 0.3;
  printf("%.40f\n", (x*x*x));

  long double y = 0.3;
  printf("%.40f\n", (double)(y*y*y));
}

这使：

0.0269999999999999996946886682280819513835
0.0269999999999999962252417162744677625597

这种情况很奇怪，因为使用更多位数的计算会产生最差的结果。这是因为无论如何初始数字 0.3 都是用几位数字近似的，因此我们从一个相对“大”的错误开始。在这种特殊情况下，发生的情况是，用少量数字计算会产生另一个“大”错误，但符号相反……因此补偿了最初的错误。相反，具有更多数字的计算会产生第二个较小的错误，但第一个仍然存在。

Answer 5

阅读 Goldberg 的“What Every Computer Scientist Should Know About Floating-Point Arithmetic”（这是 Oracle 的再版）。明白吧。浮点数不是微积分的实数。抱歉，没有可用的 TL;DR 版本。