0

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
crevell
  • 449
  • 6
  • 19
  • I'd consider the problem more to be along [these lines](https://stackoverflow.com/q/33319357), @HighPerformanceMark. And that the answer is "change `inputvariable` from `real*8` to `real`" or "change `0.1` in the source code to `0.1d0`" (both with the caveat that `real*8` means `double precision`). – francescalus May 18 '17 at 16:45
  • Thanks for the input. I worked out that the comparison was to 0.1 as `real*4` not `real*8`. What I'm wondering now is whether there is any way of testing for equality between 0.1 in `real*4` for and 0.1 in `real*8` form? It seems that there should be some way of testing that they are equal without them being the same precision. – crevell May 18 '17 at 17:04
  • 1
    Upon further digging, it seems that in general it is not advisable to test for floating point equality. I'll find a way to rewrite the code with `.LE.` instead. – crevell May 18 '17 at 17:15
  • 1
    The issue is they are not the same number when expressed in real*4 and real*8, or 0.1 when expressed as a real*8 does not equal 0.1d0. 0.1 in real*4 is 0.1000000 while in real*8 is 0.100000000000000. Note that the real*8 has the extra zeros, when the real*4 is promoted to real*8 to do the comparison the computer will fill in the missing zero's with junk, on my machine 0.1 in real*8 is 0.1000000001401161. try this by writing out a real*4 with a format specifier for double precision. – Rob May 18 '17 at 18:36
  • Very interesting that the change from `real*4` to `real*8` fills the remaining space with junk. Thanks. – crevell May 20 '17 at 08:47
  • 1
    Actually, when converting from `REAL*4` to `REAL*8`, it will fill the remaining space with zeros, not junk. However, 0.1 cannot be represented exactly in binary, so they are NOT zeros in the `REAL*8` version. Try with `0.5` or `0.25` or `0.125`, which can be exactly represented in binary, and you will see that they `.EQ.` does work for them. – Jack Jun 04 '17 at 15:56

1 Answers1

2

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.

Silasvb
  • 312
  • 1
  • 16