c - 如何在 32 位和 64 位模式下获得相同的 double 运算行为？

Question

我有一个要转换为 64 位的库。但是，我无法在 64 位模式下获得精确到位的结果，因此我的测试失败了。

我将问题简化为一个简单的测试用例:

#include <stdio.h>

int main(void) {
    printf("%d bits: ", sizeof(void*) * 8);
    volatile double d = 10.870191700000001;
    volatile double x = 0.10090000000000002;
    d += x * 30.07;
    printf("%0.15f\n", d);
}

为避免编译器差异，我使用相同的编译器和交叉编译。在这种情况下，我在 Windows 7 上的 Core i5 CPU 中使用 TDM-GCC 64 位 5.1.0。这是我的命令行:

gcc double_test.c -o double_test.exe -m32 -O0 && double_test.exe && gcc double_test.c -o double_test.exe -m64 -O0 && double_test.exe

输出是:

32 bits: 13.904254700000001
64 bits: 13.904254700000003

在这种情况下，错误很小，但在我的完整测试用例中，错误加起来足以使我的输出加倍。

如何获得与 32 位输出匹配的位精确运算？

我最接近相关的东西是使用 -ffloat-store，但在这个片段中，它获得了 32 位执行，就像 64 位执行一样，而我需要的恰恰相反。但是，这对我的图书馆没有任何明显影响。我还测试了 -fexcess-precision=standard 和 -mfp-math 选项无济于事。

Answer 1

既然你说你需要更精确的 ...01 结果，以及确定性，不幸的是你不能只使用 -msse2 -mfpmath=sse你的 32 位版本。 future 寻求确定性的读者应该使用它。

---

您可以使用 -mfpmath=387 要求 gcc 在 64 位模式下使用慢速/过时的 x87 数学，这不是默认设置。调用约定在 xmm 寄存器中传递/返回 FP args，因此这比 32 位模式更糟糕，有时需要额外的存储/重新加载。

peter@volta:/tmp$ gcc -m64 -mfpmath=387 -O3 fp-prec.c -o fp-64-387
peter@volta:/tmp$ ./fp-64-387 
64 bits: 13.904254700000001

我不确定当自动矢量化成为可能时，gcc 是否严格限制自己为 x87。如果是这样，您就错过了性能。

---

顺便说一句，在您的示例中，...01 是在 x*30.07< 的 80 位临时值中保持 extra 精度的结果 在将其添加到 d 之前。 (d 是 volatile，但是 d += stuff 仍然等同于 d = d + stuff 所以 x*30.07 不会首先四舍五入为 64 位 double)。

您可以使用long double，例如d += x * (long double)30.07 在那里强制使用 80 位临时值。 long double 在 x86-64 System V ABI 中是 80 位Linux/OS X/*BSD/etc，但在 x64 Windows 上它与 64 位 double 相同。所以这可能不适合您。

---

在这种情况下，您可以使用 FMA 获得相同的结果，它在执行加法之前保持乘法的无限精度。这在没有 FMA 支持的硬件上速度很慢，但是 fma(d, 30.07, x) 会可靠地给出您想要的结果。

如果您需要它，请在需要精度的地方使用它。

如果您在启用 FMA 的情况下进行编译，它可以内联到 FMA 指令。 (例如 -march=native 在我的 Skylake CPU 上)

即使不使用 fma() math.h 函数，gcc 也会在优化时将 mul+add 表达式收缩到 FMA 中。 (与 Clang 不同，我认为在没有 -ffast-math 的情况下，Clang 默认情况下不会执行 FP_CONTRACT)。请注意，我没有使用-march=387

 # your original source code, using an FMA instruction (native=skylake in my case)
peter@volta:/tmp$ gcc -m64 -march=native -O3 fp-prec.c -o fp-64-native
peter@volta:/tmp$ ./fp-64-native 
64 bits: 13.904254700000001

main 的相关部分是:

 57e:   c5 fb 10 44 24 08       vmovsd xmm0,QWORD PTR [rsp+0x8] # load x
 584:   c5 fb 10 0c 24          vmovsd xmm1,QWORD PTR [rsp]     # load d
 589:   c4 e2 f1 99 05 d6 01 00 00      vfmadd132sd xmm0,xmm1,QWORD PTR [rip+0x1d6]        # the 30.07 constant
 592:   c5 fb 11 04 24          vmovsd QWORD PTR [rsp],xmm0     # store d
 597:   c5 fb 10 04 24          vmovsd xmm0,QWORD PTR [rsp]     # reload d
 59c:   e8 8f ff ff ff          call   530 <printf@plt>

---

FP 确定性一般很难。

另见 https://randomascii.wordpress.com/2013/07/16/floating-point-determinism/和 https://randomascii.wordpress.com/2012/03/21/intermediate-floating-point-precision/