为什么两个相似的浮点计算会给出两个不同的结果？

Question

下面的代码通过使用特征向量作为容器或简单的C数组来实现相同的计算。它产生一个封闭的但不是位到位等效的结果。

最后的数学运算是x*alpha+y*beta。

#include <Eigen/Eigen>

int main()
{
  Eigen::VectorXd x(2);
  double* y = new double[2];
  long long int a = 4603016991731078785;
  double ga = *(double*)(&a);
  long long int b = -4617595986472363966;
  double gb = *(double*)(&b);
  long long int x0 = 451;
  long long int x1 = -9223372036854775100;
  x[0] = *(double*)(&x0);
  y[0] = *(double*)(&x0);
  x[1] = *(double*)(&x1);
  y[1] = *(double*)(&x1);
  double r = ga*x[0] + gb*x[1];
  double s = ga*y[0] + gb*y[1];
}

Answer 1

这可能是因为一个计算完全在FPU（浮点单元）内完成，精度为80位，而另一个计算使用了部分64位精度（相当于double的大小）。这也可以在不使用特征值的情况下进行演示。看下面这个程序：

int main()
{
  // Load ga, gb, y[0], y[1] as in original program
  double* y = new double[2];
  long long int a = 4603016991731078785;
  double ga = *(double*)(&a);
  long long int b = -4617595986472363966;
  double gb = *(double*)(&b);
  long long int x0 = 451;
  long long int x1 = -9223372036854775100;
  y[0] = *(double*)(&x0);
  y[1] = *(double*)(&x1);

  // Compute s as in original program
  double s = ga*y[0] + gb*y[1];

  // Same computation, but in steps
  double r1 = ga*y[0];
  double r2 = gb*y[1];
  double r = r1+r2;
}

如果你在没有优化的情况下编译这个，你会看到r和s有不同的值（至少，我在我的机器上看到了这一点）。看汇编代码，在第一次计算中，在FPU中加载ga、Y[0]、gb和Y[1]的值，然后做计算ga*Y[0]+gb*Y[1]，然后将结果存入内存。FPU用80位执行所有计算，但当结果存储在内存中时，数字会四舍五入，以便它适合于双变量的64位。

第二个计算的进行方式不同。首先，ga和Y[0]被加载到FPU中，相乘，然后四舍五入为64位数字并存储在存储器中。然后，gb和Y[1]加载到FPU中，相乘，然后舍入为64位数字并存储在存储器中。最后，r1和r2被加载到FPU中，相加，四舍五入为64位数字，并存储在存储器中。这一次，计算机轮询中间结果，这就导致了差异。

现在，我不太确定的地方来了（如果这是您的问题，我很抱歉）：这与原始程序有什么关系，其中x是一个特征容器？这里的计算是这样进行的：调用一个来自特征的函数来得到x[0]，然后ga和该函数的结果被加载到FPU中，相乘，并存储在一个临时存储器位置（64位，因此四舍五入）。然后将gb和X[1]加载到FPU中，相乘，与存储在临时存储位置中的中间结果相加，最后存储在X中。所以在原程序中r的计算中，GA*X[0]的结果被舍入到64位。这样做的原因可能是浮点堆栈没有跨函数调用保留。