ifort与gfortran计算acos(x)结果不同的原因及对齐方法

我用gfortran和ifort编译以下Fortran程序：

! acos.f90
function real8_to_int8(real8) result(int8)
    real(8), intent (in) :: real8
    integer(8) :: int8
    int8 = transfer(real8, int8)
end function

function int8_to_real8(int8) result(real8)
    integer(8), intent (in) :: int8
    real(8) :: real8
    real8 = transfer(int8, real8)
end function

program main
    integer(8) :: real8_to_int8
    real(8) :: int8_to_real8
    real :: x, acos_x

    x = int8_to_real8(4605852290655121993)

    print *, " size of x: ", sizeof(x)
    write (*, '(A, F65.60)') ' x: ', x
    print *, ""

    acos_x = acos(x)

    write (*, '(A, F65.60)') ' acos(x): ', acos(x)
    print *, "bits:     ", real8_to_int8(acos_x)
    print *, ""
    print *, ""
end program

为精确对比结果，我打印了变量的位表示，发现二者输出不同：

ifort编译运行结果：

$ ifort -real-size 64 -o acos_ifort acos.f90; ./acos_ifort

  size of x:                      8
 x:    0.852326110783516388558211929193930700421000000000000000000000
 
 acos(x):    0.550379481046229246388179490168113261461000000000000000000000
 bits:        4603132597196780746

gfortran编译运行结果：

$ gfortran -fdefault-integer-8 -fdefault-real-8 -o acos_gfortran acos.f90; ./acos_gfortran

  size of x:                     8
 x:    0.852326110783516388558211929193930700421333312988281250000000
 
 acos(x):    0.550379481046229357410481952683767303824424743652343750000000
 bits:       4603132597196780747

直接对比关键结果：

ifort:
acos(x): 0.550379481046229246388179490168113261461000000000000000000000
bits: 4603132597196780746

gfortran:
acos(x): 0.550379481046229357410481952683767303824424743652343750000000
bits: 4603132597196780747

请问：acos(x)的结果差异是否属于正常情况？或者如何修改gfortran的编译选项，使其计算结果与ifort一致？

回答 #

1. 结果差异属于正常情况 #

这种差异是浮点数计算实现细节不同导致的，完全正常。

不同编译器的数学库（ifort依赖Intel MKL，gfortran使用GNU libm）对acos这类超越函数的实现算法、精度优化策略存在差异，最终输出的结果会落在IEEE双精度浮点数的允许误差范围内（两者差异约为1e-16，远小于双精度的机器epsilon≈2.2e-16）。且两个结果都是acos(x)的正确舍入结果——只是分别舍入到了相邻的两个双精度值。

2. 让gfortran结果与ifort一致的方法 #

如果必须强制结果一致，可以尝试以下方式：

• 链接Intel MKL库替代GNU libm：
  编译时指定链接Intel MKL，让gfortran调用与ifort相同的数学函数实现。示例编译命令：

  gfortran -fdefault-integer-8 -fdefault-real-8 -o acos_gfortran acos.f90 -mkl=parallel
  

  注意需先配置好MKL的环境变量（如source Intel提供的环境脚本）。

• 调整gfortran的数学库精度选项：
  GNU libm提供部分精度控制选项，比如尝试-ffast-math（但该选项会牺牲严格IEEE兼容性，可能引发其他问题），或搭配-frounding-math -fsignaling-nans调整舍入模式，但这种方式无法保证底层算法完全匹配，不一定能实现结果完全一致。

• 手动指定舍入模式：
  在代码中通过Fortran的IEEE_SET_ROUNDING_MODE等过程强制统一舍入模式（比如默认的舍入到最近偶数），但这也无法消除算法层面的差异。

内容的提问来源于stack exchange，提问作者Rahn