Comparing floating point numbers (part 2)

From: Paul Van Delst (paul.vandelst_at_noaa.gov)
Date: 08/31/04 

Date: Tue, 31 Aug 2004 10:52:42 -0400

Hello,

Based on a question I asked the other day here, I put together the following code
(I've only shown the "double" precision code for brevity):

MODULE Compare_Float_Numbers
   USE Type_Kinds ! Contains "Single" and "Double" kind definitions
   IMPLICIT NONE
   PRIVATE
   PUBLIC :: OPERATOR (.EqualTo.)
   PUBLIC :: OPERATOR (.GreaterThan.)
   PUBLIC :: OPERATOR (.LessThan.)
   INTERFACE OPERATOR (.EqualTo.)
     MODULE PROCEDURE Is_Equal_To_Single
     MODULE PROCEDURE Is_Equal_To_Double
   END INTERFACE OPERATOR (.EqualTo.)
   INTERFACE OPERATOR (.GreaterThan.)
     MODULE PROCEDURE Is_Greater_Than_Single
     MODULE PROCEDURE Is_Greater_Than_Double
   END INTERFACE OPERATOR (.GreaterThan.)
   INTERFACE OPERATOR (.LessThan.)
     MODULE PROCEDURE Is_Less_Than_Single
     MODULE PROCEDURE Is_Less_Than_Double
   END INTERFACE OPERATOR (.LessThan.)

CONTAINS

   ELEMENTAL FUNCTION Is_Equal_To_Double( x, y ) RESULT( Equal_To )
     REAL( Double ), INTENT( IN ) :: x, y
     LOGICAL :: Equal_To
     Equal_To = ABS( x - y ) < SPACING( MAX(ABS(x),ABS(y)) )
   END FUNCTION Is_Equal_To_Double

   ELEMENTAL FUNCTION Is_Greater_Than_Double( x, y ) RESULT ( Greater_Than )
     REAL( Double ), INTENT( IN ) :: x, y
     LOGICAL :: Greater_Than
     IF ( (x - y) > SPACING( MAX( ABS(x), ABS(y) ) ) ) THEN
       Greater_Than = .TRUE.
     ELSE
       Greater_Than = .FALSE.
     END IF
   END FUNCTION Is_Greater_Than_Double

   ELEMENTAL FUNCTION Is_Less_Than_Double( x, y ) RESULT ( Less_Than )
     REAL( Double ), INTENT( IN ) :: x, y
     LOGICAL :: Less_Than
     IF ( (y - x) > SPACING( MAX( ABS(x), ABS(y) ) ) ) THEN
       Less_Than = .TRUE.
     ELSE
       Less_Than = .FALSE.
     END IF
   END FUNCTION Is_Less_Than_Double
END MODULE Compare_Float_Numbers

and I can use the functions like so

   REAL :: x, y
   ....

   IF ( x .GreaterThan. y ) THEN
    ....
   END IF

Now, I'm testing this code by setting (for the double precision case),
   x = ..some number..
   y1 = NEAREST( x, 1.0_Double )
   y2 = y1 - SPACING( x )

Some of the numbers generated look like:

      Double test. x = 1.23456789012345603894E-16
                   y1 = 1.23456789012345628545E-16 : NEAREST( x, 1.0 )
                   y2 = 1.23456789012345603894E-16 : y1 - SPACING( x )

      Double test. x = 1.23456789012345594103E-01
                   y1 = 1.23456789012345607981E-01 : NEAREST( x, 1.0 )
                   y2 = 1.23456789012345594103E-01 : y1 - SPACING( x )

      Double test. x = 1.23456789012345600000E+16
                   y1 = 1.23456789012345620000E+16 : NEAREST( x, 1.0 )
                   y2 = 1.23456789012345600000E+16 : y1 - SPACING( x )

etc.. where to the resolution printed out, x < y1 and x == y2

The results I get for the .EqualTo. operator are as expected
   x .EqualTo. y1 = F
   x .EqualTo. y2 = T

and for the other operators I always get:
   x .GreaterThan. y1 = F
   x .GreaterThan. y2 = F
   x .LessThan. y1 = F
   x .LessThan. y2 = F

I don't understand the x .LessThan. y1 result. The .EqualTo. tells me they're not equal,
and the way I constructed the numbers above gives me a y1 that is *always* greater than x, so
how come the x .LessThan. y1 is always false? Is it a problem with how I've defined the
functions, or how I've constructed the test numbers? Or am I (again) misunderstanding how
floating point numbers "work"?

Any insights, suggestions are welcome.

cheers,

paulv