如何在32位和64位模式下获得双精度操作的相同行为？

我有一个库,我正在转换为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);
}

为避免编译器差异,我使用相同的编译器和交叉编译.在这种情况下,我在Core i5 cpu上使用Windows 7上的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,但是在这个片段中,它得到了像64位一样的32位执行,而我需要恰恰相反.但是,这对我的图书馆没有任何明显的影响.我还测试了-fexcess-precision = standard和-mfp-math选项无济于事.

解决方法

既然你说你需要更精确的… 01结果,以及确定性,你遗憾的是你不能在你的32位构建中使用-msse2 -mfpmath = sse.寻找决定论的未来读者应该使用它.

您可以使用-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.如果是这样,你会错过表现.

而BTW,在你的例子中,… 01是在将x + 30.07添加到d之前为x * 30.07保持80位临时值的额外精度的结果. (d是易失性的,但是d = stuff仍然相当于d = d的东西,因此x * 30.07不会首先四舍五入到64位双精度).

你可以使用long double,例如d = x *(long double)30.07强制80位临时存在.在x86-64 System V ABI Linux / OS X / * BSD /等中,long double是80位,但在x64 Windows上,它与64位double相同.所以这可能不适合你.

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

如果需要,请在需要精度的地方使用它.

如果在启用FMA的情况下进行编译,则可以内联到FMA指令. (例如-march =我的Skylake cpu上的原生)

即使不使用fma()math.h函数,gcc也会在优化时将mul添加表达式添加到FMA中. (与Clang不同,我认为默认情况下不使用-ffast-math进行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

主要的相关部分是：

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/

如何在32位和64位模式下获得双精度操作的相同行为？

解决方法

相关文章