CODEWITHSUNDEEP
fortran
crevell
crevell

Reputation: 460

Fortran floating point equality

I have a Fortran program that tests equality in two floating point numbers. It can be condensed to what is shown below. When this program is run with "0.1" given as a command line argument, I expect it to print "what I expected" but instead it prints "strange". I understand that this is probably due to a floating point rounding issue, but am hoping someone might be able to explain exactly how I should change inputvariable to make this code print "what I expected" with a command line argument of 0.1

program equalitytest
  character(len=3) :: arg1
  real*8           :: inputvariable
  CALL GET_COMMAND_ARGUMENT(1,arg1)
  READ(arg1,*) inputvariable
  IF (inputvariable.EQ.0.1) THEN
    PRINT*, "what I expected"
  ELSE
    PRINT*, "strange"
  ENDIF
end program equalitytest

Run as follows:

./equalitytest 0.1
strange

Upvotes: 0

Views: 1690

Answers (1)

Silasvb
Silasvb

Reputation: 332

As a general point, there should be very few reasons why one should need to compare equality with real numbers. If I ever find myself writing such code, I tend to pause and have a think about what I am trying to achieve. What real-world condition is actually a reflection of this?

The exception to the above relate to zeros, either when writing robust code which checks for and handles possible divisions by zero, or for cases searching for a convergent solution to an equation - in the latter case, this should be considered using a delta anyway.

If there really is a need for this check, why not outsource it to a standard library within the project, e.g.

module mylib
    use iso_fortran_env

    implicit none
    private

    public :: isReal4EqualReal4
    public :: isReal4EqualReal8
    public :: isReal8EqualReal4
    public :: isReal8EqualReal8

    real(real32), parameter :: delta4 = 0.001
    real(real64), parameter :: delta8 = 0.0000000001

    contains

        logical function isReal4EqualReal4(lhs, rhs) result(equal)
            real(real32), intent(in) :: lhs
            real(real32), intent(in) :: rhs
            equal = (abs(lhs - rhs) .le. delta4)
        end function isReal4EqualReal4

        logical function isReal4EqualReal8(lhs, rhs) result(equal)
            real(real32), intent(in) :: lhs
            real(real64), intent(in) :: rhs
            equal = (abs(lhs - real(rhs,4)) .le. delta4)
        end function isReal4EqualReal8

        logical function isReal8EqualReal4(lhs, rhs) result(equal)
            real(real64), intent(in) :: lhs
            real(real32), intent(in) :: rhs
            equal = isReal4EqualReal8(rhs, lhs)
        end function isReal8EqualReal4

        logical function isReal8EqualReal8(lhs, rhs) result(equal)
            real(real64), intent(in) :: lhs
            real(real64), intent(in) :: rhs
            equal = (dabs(lhs - rhs) .le. delta8)
        end function isReal8EqualReal8

end module mylib

EDIT: Forgot to add that one of the benefits of the above is the compiler will warn me if I'm attempting to compare two real numbers of different types while using the wrong interface

EDIT: Updated to use portable real number definitions.

Upvotes: 2

Related Questions