Geek Answers Handbook

Floating Point Equality Testing

I use gfortran in MinGW under Windows 7 (32 bit) to compile Fortran code. Here is the minimum code contained in the testequal.f file:

  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

output:

 1.0000000000000000 1.0000000000000000 F F

But I expect the two gates to be like T (true). What is the problem?

+7
fortran mingw
source share
3 answers

With rare exceptions, do not compare floating point numbers for exact equality. Finite precision floating point arithmetic rules are not the same as real number arithmetic rules. Compare numbers with tolerance e.g.

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

The correct comparison for realities should be agnostic. This is a good way to compare:

 if (abs(a1-a2) <= epsilon(a1)) print*, 'a1=a2'
+1
source share

The real problem is that 0.3 and 0.7 cannot be exactly expressed in binary format.
0.3 => 0.010011001100110011 ....... 0.7 => 0.0101100110011001100 .......

When they are saved, if both are rounded or both are rounded down, adding two numbers will not return to 1.000000000 .....

This is a common mistake in modeling programming. Try to have step sizes that are natural for a computer:

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

The negative powers of the two can be represented exactly on the computer. So it really works:

 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
0
source share

More articles:


All Articles