Windows 下 Visual Studio 2012 / Intel 编译器的 C++ 双精度失败答案

【问题标题】：C++ double precision failure with Visual Studio 2012 / Intel Compiler under WindowsWindows 下 Visual Studio 2012 / Intel 编译器的 C++ 双精度失败
【发布时间】：2014-05-10 05:50:45
【问题描述】：

我在使用 Visual Studio 和英特尔编译器的 Windows 下遇到双精度问题。

#include <stdio.h>
#include <math.h>
int main()
{
  double result = 42;
  int rounds = 14;
  for(int a=1; a<=rounds; a++)
  {
       result = sqrt(sqrt(result));
       printf("Round %2.i is %.17g\n", a, result);
   }
   for(int a=1; a<=rounds; a++)
   {
       result = pow(result, 4);
       printf("Round %2.i is %.17g\n", a, result);
   }
  return 0;
}

使用这些代码，我在 Linux 下使用这些编译器设置生成：

linux下的Intel编译器：

icc -O2 -fp-model double test.cpp -o test

linux下的GCC编译器：

g++ -O2 test.cpp -o test

和Windows下的Matlab和Scilab的正确结果：

Round  1 is 2.5457298950218306
Round  2 is 1.2631446315995756
Round  3 is 1.0601401197016873
Round  4 is 1.0147073765340913
Round  5 is 1.0036567375971333
Round  6 is 1.0009129334669549
Round  7 is 1.0002281552726189
Round  8 is 1.0000570339386641
Round  9 is 1.0000142581797196
Round 10 is 1.0000035645258711
Round 11 is 1.0000008911302767
Round 12 is 1.0000002227824947
Round 13 is 1.000000055695619
Round 14 is 1.0000000139239045
Round  1 is 1.000000055695619
Round  2 is 1.0000002227824947
Round  3 is 1.0000008911302765
Round  4 is 1.0000035645258705
Round  5 is 1.0000142581797171
Round  6 is 1.0000570339386543
Round  7 is 1.0002281552725802
Round  8 is 1.0009129334668005
Round  9 is 1.0036567375965142
Round 10 is 1.0147073765315875
Round 11 is 1.0601401196912235
Round 12 is 1.263144631549705
Round 13 is 2.545729894619797
Round 14 is 41.999999973468661

在 linux 下使用这个 Intel 编译器选项

icc -O2 test.cpp -o test

linux下的这个gcc编译器选项

g++ -O2 -ffast-math test.cpp -o test

和所有 Visual Studio 2012 编译器设置我得到这个错误的结果：

Round  1 is 2.5457298950218306
Round  2 is 1.2631446315995756
Round  3 is 1.0601401197016873
Round  4 is 1.0147073765340913
Round  5 is 1.0036567375971333
Round  6 is 1.0009129334669549
Round  7 is 1.0002281552726189
Round  8 is 1.0000570339386641
Round  9 is 1.0000142581797196
Round 10 is 1.0000035645258711
Round 11 is 1.0000008911302767
Round 12 is 1.0000002227824947
Round 13 is 1.000000055695619
Round 14 is 1.0000000139239045
Round  1 is 1.000000055695619
Round  2 is 1.0000002227824947
Round  3 is 1.0000008911302767
Round  4 is 1.0000035645258711
Round  5 is 1.0000142581797198
Round  6 is 1.0000570339386647
Round  7 is 1.0002281552726222
Round  8 is 1.0009129334669684
Round  9 is 1.0036567375971872
Round 10 is 1.0147073765343093
Round 11 is 1.0601401197025984
Round 12 is 1.2631446316039172
Round 13 is 2.5457298950568319
Round 14 is 42.000000002309839

在 Windows 下，以下设置被忽略，我总是得到错误的结果。

windows下的Intel编译器选项：

icl /O2 /fp:double test.cpp -o test.exe

Windows 下的 Visual Studio 编译器选项：

/fp:precise
/fp:strict
/fp:fast

产生了所有相同的“错误数字”......

为什么会出现这个错误，我该如何解决我在所有平台上获得相同工作双精度的问题？

非常感谢和问候

娜塔莉

编辑：

为了进行基准测试，我循环了代码：

#include <stdio.h>
#include <math.h>
int main()
{
  double x = 1;
  double result = 0;
  int rounds = 12;
  while (x!=result)
  {
  x++;
  result=x;
  for(int a=1; a<=rounds; a++)
  {
       result = sqrt(sqrt(result));
   }
   for(int a=1; a<=rounds; a++)
   {
       result = pow(result, 4);
   }
   }
   printf("The next possible with %2.i rounds is %.0f\n", rounds, result);
  return 0;
}

使用工作的 linux 编译器编译：英特尔编译器：

icc -O2 -fp-model double speedtest3.cpp -o speedtest3

结果：

The next possible with 12 rounds is 3671078
real    0m1.950s
user    0m1.947s
sys     0m0.003s

GCC 编译器：

g++ -O2 speedtest3.cpp -o speedtest3

结果：

The next possible with 12 rounds is 3671078

real    0m3.445s
user    0m3.442s
sys     0m0.004s

带有 pow 的代码：

#include <stdio.h>
#include <math.h>
int main()
{
  double x = 1;
  double result = 0;
  int rounds = 12;
  while (x!=result)
  {
  x++;
  result=x;
  for(int a=1; a<=rounds; a++)
  {
       result = pow(result, 0.25);
   }
   for(int a=1; a<=rounds; a++)
   {
       result = pow(result, 4);
   }
   }
   printf("The next possible with %2.i rounds is %.0f\n", rounds, result);
  return 0;
}

使用相同的选项编译：结果英特尔编译器：

The next possible with 12 rounds is 3671078

real    0m2.887s
user    0m2.885s
sys     0m0.004s

结果 GCC 编译器：

The next possible with 12 rounds is 3671078

real    0m5.905s
user    0m5.905s
sys     0m0.008s

结果：

sqrt 比 pow 快近 2 倍 intel 编译器比 gcc 编译器快近 2 倍

【问题讨论】：

至少对我来说，不清楚您尝试过哪些平台/编译器，哪些失败了，哪些是正确的。
42.000000002309839 错误和 41.999999973468661 正确的依据是什么？例如，第一个比第二个稍微接近 42，尽管考虑到误差的指数，它们都可能是合理的。顺便说一句，在任何 C 实现中，result*=result;result*=result 甚至result=result*result*result*result 都可能比result=pow(result,4) 更快、更准确，而pow(result, 0.25) 很可能比sqrt(sqrt(result)) 更好。
感谢您的回答。我希望 Scilab 和 Matlab 的结果在“双精度”（41.999999973468661）下是正确的。
据我了解，42.000000002309839 必须更接近正确的结果。这是长双精度吗？我很快就会对不同的方式进行性能测试。平台：Linux：Gentoo Kernel 3.13.6 强化 x86_64、gcc 4.7.3、icc 14.0.2 Windows：Windows 7 Pro 64bit、Visual Studio 2012 最新更新、icc 14.0.3.202
使用result=result*result*result*result 我得到一个全新的不同号码：42.000000637677026 但pow(result, 0.25) 工作正常。结果是一样的。

标签： c++ visual-studio gcc compiler-construction double-precision

【解决方案1】：

他们都是正确的。

在处理任何精度的浮点数学时，您不能要求结果是任何精确到最后一位数字的数字。期望准确地得到 41.999999973468661 与期望准确地得到 42.000000000000000 一样被误导。

下面是一个示例，说明如何获得不同的结果：C++ 中的双精度浮点运算是（通常取决于编译器和平台）64 bits。但是，英特尔 CPU 的浮点硬件在内部使用 80 bits 计算数字。在某些时候，这些值会从 FPU 中提取并修整为 64 位。 “在某些时候”的含义取决于编译器及其优化器认为什么是合适或最佳的，因此，不同的编译器/设置将导致计算的不同部分四舍五入，从而影响结果的确切数字。

即使在不同的计算机上使用完全相同的二进制程序也可能合法地改变结果，并且仍然是“正确的”。例如，过去某些英特尔处理器根本没有硬件浮点单元（很久以前，最后一个是 486SX）。在这些情况下，您的程序将在一个单独的纯软件库上运行计算，该库在 80 位精度下什么都不做。

附录：

我想正确看待事物。 Wolfram Alpha 说纽约和旧金山之间的距离约为 2577 英里。在那个尺度上，你如此担心两个数字之间的差异就像说你自己对距离的估计“错了 1/10^th 英寸”^{(*)支持>。如果这种差异对您的问题很重要，那么浮点运算就不是适合这项工作的工具。}

^* 相对误差为 (42.000000002309839 - 41.999999973468661)/41.999999973468661 = 6.8669471471949833966026905617771e-10。 2577 英里 = 163278720 英寸。该测量值的等效差异为 (163278720 英寸 * 6.8669471471949833966026905617771e-10) = 0.1121226...英寸。所有计算都在 Windows 7 上的 Windows Calculator 中完成，它使用定点算法。

【讨论】：

谢谢，有趣的一点。有没有办法使用 80 位硬件并将计算的结果减少到 64 位？在工作中，数学专家使用 matlab，而我必须使用 scilab。但是我可以使用 c++ 比 matlab 和 scilab 更快地编写代码。但我必须确保数值计算结果与 matlab/scilab 中的相同。在家里我可以使用 linux，但在工作中我必须使用 Windows。
@Natalie，你正在打一场失败的战斗。 FPU 的内部精度只是影响精确结果的众多因素之一。这是浮点运算的一个基本限制，您应该从不比较两个值是否相等，只有它们的差异是否超过了一个实质性的容差（对您的特定问题有实质性);不在 C++ 中，甚至在 MathLab 中也没有。每项技术都是一种妥协：在 FP 数学中，我们用无限的精度换取可预测的速度。由于 MathLab 本身使用 FP，因此无论您使用什么工具，都会受到这种妥协的约束。

【解决方案2】：

看看 pow(double,double) 的一种实现方式，here。来自不同供应商的不同编译器必然会包含不同的实现，从而产生差异。

正如建议的那样，更好的选择是使用

double h = x*x;
x = h*h;

你所说的“不正确”更准确，比较一下：

delta is 2.3098394308362913e-09   42.00000000230983 - 42
delta is 2.6531338903623691e-08   42 - 41.999999973468661

稍后你想证明或展示什么？重复执行一些数值计算通常比直接计算结果更糟糕；因此 pow( 42, 1.0/(1024*1024*64)) 理论上比在循环中取第四个根更可取，相反的过程也是如此： pow( result, 1024*1024*64 )，这可能使 pow 能够使用最少的乘法。

无论如何，42 的无数次方根不能表示为机器数，因此它的无数次方必然会偏离 42。您可以使用初等数值数学来估计误差。

【讨论】：

您的代码 double h = x*x; 可以正常工作并生成更准确的 42.000000002309839。但在工作中，数学专家使用 matlab，我必须使用 scilab。但是我可以用 c++ 比 matlab 和 scilab 更快地编写代码。但我必须确保数值计算结果与 matlab/scilab 中的相同。在家里我可以使用 linux，但在工作中我必须使用 Windows。