CODEWITHSUNDEEP
fortranmingw
Hao Wang
Hao Wang

Reputation: 519

Testing floating point numbers equality

I am using gfortran in MinGW under Windows 7 (32bit) to compile Fortran code. Here is the minimal code contained in the file 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

Compiled with

gfortran testequal.f -std=legacy

the output is:

1.0000000000000000
1.0000000000000000
F
F

But I expect the two booleans to be both T (true). What is the problem here?

Upvotes: 8

Views: 9546

Answers (3)

Kostas
Kostas

Reputation: 81

A proper comparison for reals should be agnostic of kind. This is a good way to compare:

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

Upvotes: 4

Jack
Jack

Reputation: 6158

The real issue is that 0.3 and 0.7 cannot be expressed exactly in binary.
0.3 => 0.010011001100110011....... 0.7 => 0.101100110011001100.......

When they are stored, if both are rounded up or both are rounded down, the addition of the two numbers will not get back to 1.000000000.....

This is a common error in programming simulations. Try to have step sizes that are natural to the computer:

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

Negative powers of two can be represented exactly in the computer. So this will actually work:

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

Upvotes: 1

M. S. B.
M. S. B.

Reputation: 29391

With rare exceptions, don't compare floating point numbers for exact equality. The rules of finite-precision floating point arithmetic are not the same are the rules of real number arithmetic. Compare numbers with a tolerance, e.g., 

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

Upvotes: 11

Related Questions