测试浮点数相等

Posted 2021-04-16

我在Windows 7（32位）下的MinGW中使用gfortran来编译Fortran代码。这是testequal.f文件中包含的最小代码：

      program testequal
      real*8 a1, a2

      a1 = 0.3d0
      a2 = 0.7d0

      write(*,*) 1.d0
      write(*,*) a1+a2
      write(*,*) a1+a2.eq.1.0
      write(*,*) a1+a2.eq.1.d0
      end

编译

gfortran testequal.f -std=legacy

输出是：

1.0000000000000000
1.0000000000000000
F
F

但我希望两个布尔都是T（真）。这里有什么问题？

答案

除了极少数例外，不要将浮点数进行比较以确保完全相等。有限精度浮点运算的规则与实数运算的规则不一样。比较数字与容差，例如，

sum = a1 + a2
if ( abs (sum - 1.0) < 1.0D-5 ) ...

另一答案

对于实际的适当比较应该是不可知的。这是比较的好方法：

if (abs(a1-a2) <= epsilon(a1)) print*, 'a1=a2'

另一答案

真正的问题是0.3和0.7不能用二进制表示。 0.3 => 0.010011001100110011 ....... 0.7 => 0.101100110011001100 .......

当它们被存储时，如果两者都被向上舍入或两者都被向下舍入，则两个数字的相加将不会回到1.000000000 .....

这是编程模拟中的常见错误。尝试使用计算机自然的步长：

real*8 :: dt=2.0d0**(-5)

两个的负幂可以在计算机中准确表示。所以这实际上会起作用：

program negative_powers
  real*8 :: dt = 2.0d0**(-8)
  real*8 :: t = 0.0d0

  do while (t .ne. 500.0d0)
     print *, t
     t = t + dt
  end do
  print *, t
end program negative_powers