1

以下代码通过仅使用特征向量作为容器或简单的 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];
}

为什么会这样？

使用 MSVC 和 gcc（64 位操作系统）时结果不同。

4

1 回答 1

1

这可能是因为一个计算完全在 FPU（浮点单元）内完成，精度为 80 位，而另一个计算部分使用 64 位精度（双精度的大小）。这也可以在不使用 Eigen 的情况下证明。看下面的程序：

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 有不同的值（至少，我在我的机器上看到了）。看汇编代码，在第一次计算中，将ga、y[0]、gb和y[1]的值加载到FPU中，那么计算ga * y[0] + gb * y[1]为完成，然后将结果存储在内存中。FPU 使用 80 位进行所有计算，但是当结果存储在内存中时，数字会四舍五入，以使其适合双精度变量的 64 位。

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

对于此计算，舍入具有相当大的影响，因为您正在使用非正规数。

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

于 2013-02-04T12:39:33.570 回答