Topic: No Intrinsic Procedure for real equals real?

Hi Jeff,

I have a general programming question regarding the precision of real variables. Particularly when testing two reals for equality (x == y).  Due to the lack of precision, there are situations when x should equal y, but due to the error in precision they are not equal (1.0000000000000000 /= 0.9999999999999999).  As far as I understand, there is not a FORTRAN intrinsic procedure for the equality of two real numbers other than x == y.  Therefore, I've written a procedure to test for the equality of two real numbers that allows for the lack of precision in variables x and y.  The routine is copied below.  My question is as follows, (1) Is this the best way to test for equality?, (2) Is there a better approach that works well?, (3) Why does double precision x+factor == y1 (.true.) in my example?  I welcome any thoughts on or improvements in this approach.  Finally, Jeff, if you think my question is not appropriate for this forum, please remove this message.

Frank

Equality Program:

        PROGRAM Test_EQ_FP8
        !---------------------------------------------------------------
        ! NAME:
        !       Test_EQ_FP8
        ! PURPOSE:
        !       Program to test the routines in the EQ_FP8(x,y) routine.
        ! MODULES:
        !       Type_Kinds: Definitions for kinds of variable types.
        !       EQ_FP8:     Performs equality comparison of real numbers.
        !---------------------------------------------------------------

        IMPLICIT NONE
        !-------------------------------
        !--- Floating Point Precision
        !-------------------------------
        !    REAL( KIND )
        !===============================
        !    Single (4  bytes)
        !    Double (8  bytes)
        !    Quad   (16 bytes)
        !-------------------------------
        !--- Six decimal place precision
        INTEGER, PARAMETER :: &
        Single = SELECTED_REAL_KIND(6)

        !--- Fourteen decimal place precision
        INTEGER, PARAMETER :: &
        Double = SELECTED_REAL_KIND(2*precision(1.0_Single))

        !--- Sixteen decimal place precision
        INTEGER, PARAMETER :: &
        Quadr = SELECTED_REAL_KIND(2*precision(1.0_Double))
        INTEGER, PARAMETER :: &
        Quad = (((1 + SIGN(1,Quadr)) / 2) * Quadr) &
        + (((1 - SIGN(1,Quadr)) / 2) * Double)

        INTEGER, PARAMETER :: &
        Q = Quad, D = Double, S = Single

        ! ----------
        ! Parameters
        CHARACTER( * ), PARAMETER :: ROUTINE = 'Test_EQ_FP8'

        ! -- The test numbers
        INTEGER, PARAMETER :: N_NUMBERS = 7

        REAL( Single ), PARAMETER, DIMENSION( N_NUMBERS ) :: SINGLE_NUMBER = &
        (/ 1.000000000000000e-32_Single, &
        1.234567890123456e-16_Single, 1.234567890123456e-08_Single, &
        1.000000000000000e+00_Single, 1.234567890123456e+08_Single, &
        1.234567890123456e+16_Single, 1.234567890123456e+32_Single /)

        REAL( Double ), PARAMETER, DIMENSION( N_NUMBERS ) :: DOUBLE_NUMBER = &
        (/ 1.000000000000000e-32_Double, &
        1.234567890123456e-16_Double, 1.234567890123456e-08_Double, &
        1.000000000000000e+00_Double, 1.234567890123456e+08_Double, &
        1.234567890123456e+16_Double, 1.234567890123456e+32_Double /)

        REAL( Quad ), PARAMETER, DIMENSION( N_NUMBERS ) :: QUAD_NUMBER = &
        (/ 1.000000000000000e-32_Quad, &
        1.234567890123456e-16_Quad, 1.234567890123456e-08_Quad, &
        1.000000000000000e+00_Quad, 1.234567890123456e+08_Quad, &
        1.234567890123456e+16_Quad, 1.234567890123456e+32_Quad /)

        ! ---------
        ! Variables
        ! ---------
        CHARACTER(len=256) :: cend
        REAL( Single ) :: x,  y1, y2, y3, y4
        REAL( Double ) :: xd, yd1, yd2, yd3, yd4
        REAL( Quad   ) :: xq, yq1, yq2, yq3, yq4
        REAL( Double ) :: x0, y0
        INTEGER :: i

        !--- Added error
        REAL( Double ), PARAMETER :: e = 1.0e-30_Double

        !--- TEST EQUALITY OF X AND Y REALS
        !---------------------------------------------------------------
        DO i = 1, N_NUMBERS
            WRITE( *, '(/5x, "NUMBER ", i2, " of ", i2 )' ) i, N_NUMBERS

            !----------------------
            ! Single precision test
            !----------------------
            x = SINGLE_NUMBER(i)
            y1 = NEAREST( x, 1.0_Single )
            y2 = y1 - SPACING( x )
            y3 = NEAREST( x, -1.0_Single )
            y4 = y3 + SPACING( x )
            x = x + e

            !--- Test standalone single precision function routines
            WRITE( *, '( /5x, "SINGLE PRECISION: x  = ", es20.13, &
            &/5x, "y1 = ", es20.13, 2x, ":  NEAREST( x, 1.0 )", &
            &/5x, "y2 = ", es20.13, 2x, ":  y1 - SPACING( x )", &
            &/5x, "y3 = ", es20.13, 2x, ":  NEAREST( x,-1.0 )", &
            &/5x, "y4 = ", es20.13, 2x, ":  y3 + SPACING( x )" )' ) &
            x, y1, y2, y3, y4

            x0 = x; y0 = y1
            WRITE( *, '(/5x,a,2x,L6)') 'EQ_FP8(x,y1): X=Y?', EQ_FP8( x0, y0 )
            x0 = x; y0 = y2
            WRITE( *, '( 5x,a,2x,L6)') 'EQ_FP8(x,y2): X=Y?', EQ_FP8( x0, y0 )
            x0 = x; y0 = y3
            WRITE( *, '( 5x,a,2x,L6)') 'EQ_FP8(x,y3): X=Y?', EQ_FP8( x0, y0 )
            x0 = x; y0 = y4
            WRITE( *, '( 5x,a,2x,L6)') 'EQ_FP8(x,y4): X=Y?', EQ_FP8( x0, y0 )

            !----------------------
            ! Double precision test
            !----------------------
            xd = DOUBLE_NUMBER(i)
            yd1 = NEAREST( xd, 1.0_Double )
            yd2 = yd1 - SPACING( xd )
            yd3 = NEAREST( xd, -1.0_Double )
            yd4 = yd3 + SPACING( xd )
            xd = xd + e

            !--- Test stand alone double precision function routines
            WRITE( *, '( /5x, "DOUBLE PRECISION: x  = ", es20.13, &
            &/5x, "y1 = ", es20.13, 2x, ":  NEAREST( x, 1.0 )", &
            &/5x, "y2 = ", es20.13, 2x, ":  y1 - SPACING( x )", &
            &/5x, "y3 = ", es20.13, 2x, ":  NEAREST( x,-1.0 )", &
            &/5x, "y4 = ", es20.13, 2x, ":  y3 + SPACING( x )" )' ) &
            xd, yd1, yd2, yd3, yd4

            x0 = xd; y0 = yd1
            WRITE( *, '(/5x,a,2x,L6)') 'EQ_FP8(x,y1) X=Y?:', EQ_FP8( x0, y0 )
            x0 = xd; y0 = yd2
            WRITE( *, '( 5x,a,2x,L6)') 'EQ_FP8(x,y2) X=Y?:', EQ_FP8( x0, y0 )
            x0 = xd; y0 = yd3
            WRITE( *, '( 5x,a,2x,L6)') 'EQ_FP8(x,y3) X=Y?:', EQ_FP8( x0, y0 )
            x0 = xd; y0 = yd4
            WRITE( *, '( 5x,a,2x,L6)') 'EQ_FP8(x,y4) X=Y?:', EQ_FP8( x0, y0 )

            !-------------------------
            ! Quadruple precision test
            !-------------------------
            xq = QUAD_NUMBER(i)
            yq1 = NEAREST( xq, 1.0_Quad )
            yq2 = yq1 - SPACING( xq )
            yq3 = NEAREST( xq, -1.0_Quad )
            yq4 = yq3 + SPACING( xq )
            xq = xq + e

            !--- Test stand alone quadruple precision function routines
            WRITE( *, '( /5x, "QUADRULE PRECISION: x  = ", es20.13, &
            &/5x, "y1 = ", es20.13, 2x, ":  NEAREST( x, 1.0 )", &
            &/5x, "y2 = ", es20.13, 2x, ":  y1 - SPACING( x )", &
            &/5x, "y3 = ", es20.13, 2x, ":  NEAREST( x,-1.0 )", &
            &/5x, "y4 = ", es20.13, 2x, ":  y3 + SPACING( x )" )' ) &
            xq, yq1, yq2, yq3, yq4

            x0 = xq; y0 = yq1
            WRITE( *, '(/5x,a,2x,L6)') 'EQ_FP8(x,y1) X=Y?:', EQ_FP8( x0, y0 )
            x0 = xq; y0 = yq2
            WRITE( *, '( 5x,a,2x,L6)') 'EQ_FP8(x,y2) X=Y?:', EQ_FP8( x0, y0 )
            x0 = xq; y0 = yq3
            WRITE( *, '( 5x,a,2x,L6)') 'EQ_FP8(x,y3) X=Y?:', EQ_FP8( x0, y0 )
            x0 = xq; y0 = yq4
            WRITE( *, '( 5x,a,2x,L6)') 'EQ_FP8(x,y4) X=Y?:', EQ_FP8( x0, y0 )

            WRITE( *,'(5x,a)') REPEAT("-",50)

            !----------------------
            !--- Question?
            !----------------------
            WRITE( *, '(/5x,a/)') &
                "Why is double precision number x+e=y1 true?"
        END DO

        WRITE( *, '( /10x, "Press <ENTER> to continue..." )' )
        READ( *, '(a)' ) cend

        !---------------------------------------------------------------
        CONTAINS
        !---------------------------------------------------------------

        FUNCTION EQ_FP8( x,y ) Result( DrFP )
        !***************************************************************
        !---    EQ_FP8 version 1.0
        !---
        !---    Developed by F.W. Perrella, Ph.D.  November 7, 2020
        !---    EQ_FP8™ Copyright © 2020 F.W. PERRELLA, Ph.D.
        !---            Silver Spring, Maryland U.S.A.
        !---                All Rights Reserved.
        !***************************************************************
        !
        ! NAME:
        !       EQ_FP8
        ! PURPOSE:
        !       Program to determine whether two real numbers are equal.
        ! MODULES:
        !       Type_Kinds:      Definition for double kind variable types.
        !       EQ_FP8:         Routine to perform equality comparison
        !                       of real numbers.
        !---------------------------------------------------------------
        !
        !    SOURCE listing of EQ_FP8.f90:
        !       u = DMIN1(DABS(x), DABS(y))
        !       v = DMAX1(DABS(x), DABS(y))
        !       Difx  = NEAREST( u, v )
        !       Dify  = SPACING( v )
        !       Difxy = NEAREST( v-u, v )
        !       IF (Difxy < TINY(u)) Difxy = TINY(u)
        !       IF ( Difxy < DABS(Difx + Dify - x) ) THEN
        !           DrFP = .TRUE.
        !       ELSE
        !           DrFP = .FALSE.
        !       ENDIF
        !
        !    FORTRAN 90 INTRINSIC FUNCTIONS:
        !      SPACING(x) determines the absolute distance between the
        !      argument x and the nearest adjacent number of the same type.
        !
        !      SPACING(x) = 2.D0**( EXPONENT(x) - DIGITS(x) )
        !
        !      NEAREST( x,s ) determines the processor-representable number
        !      nearest to x in the direction indicated by the sign of s.
        !
        !      NEAREST(x,y) = x + SIGN(1.D0,x)*2.D0**(EXPONENT(x)-DIGITS(x))
        !
        !    CALLING SEQUENCE:
        !       Result = EQ_FP8( x,y )
        !
        !    INPUT ARGUMENTS:
        !       x, y:       Two real numbers to compare.
        !       TYPE:       REAL( Double )
        !
        !    ATTRIBUTES:
        !       INTENT( IN )
        !
        !
        !    FUNCTION RESULT:
        !       TYPE:     Logical
        !       RESULT:   A returned logical value indicating whether the
        !                 x and y inputs are equal within the set precision.
        !                 Set to .TRUE. if the real numbers are equal
        !                 within the specified tolerance.
        !                 Set to .FALSE. if the real numbers are not equal.
        !---------------------------------------------------------------
            IMPLICIT NONE

            !--- Input variables
            REAL(Kind=8), INTENT( IN )  :: x
            REAL(Kind=8), INTENT( IN )  :: y
            LOGICAL                     :: DrFP

            !--- Double precision
            INTEGER, PARAMETER :: Doubles = SELECTED_REAL_KIND(15,300)

            !--- Local Variables
            REAL(Doubles)            :: u, v
            REAL(Doubles)            :: Difx
            REAL(Doubles)            :: Dify
            REAL(Doubles)            :: Difxy
            INTRINSIC                :: NEAREST, SPACING, DMIN1, DMAX1, DABS

            !--- Minimum and maximum of x and y
            u = DMIN1(DABS(x), DABS(y))
            v = DMAX1(DABS(x), DABS(y))

            !--- Determine Nearest, Spacing, and Nearest difference
            Difx  = NEAREST( u, v )
            Dify  = SPACING( v )
            Difxy = NEAREST( v-u, v )

            !--- Nearest difference cannot be less than Tiny
            IF (Difxy < TINY(u)) Difxy = TINY(u)

            !--- Evaluate if x is equal to y?
            IF ( Difxy < DABS(Difx + Dify - x) ) THEN
                DrFP = .TRUE.
            ELSE
                DrFP = .FALSE.
            ENDIF

        END FUNCTION EQ_FP8
        !---------------------------------------------------------------

        END PROGRAM Test_EQ_FP8
        !---------------------------------------------------------------
        !---------------------------------------------------------------

Re: No Intrinsic Procedure for real equals real?

Hi Jeff,

The equality routine FUNCTION EQ_FP8( x,y ) previously listed should be corrected to -

        !--- Evaluate if x is equal to y?
        IF (DABS(Difxy) < TINY(u)) Difxy = TINY(u)
        IF ( x == y ) THEN
            DrFP = .TRUE.
        ELSEIF ( DABS(Difxy) < 2.D0*DABS(Difx + Dify - DABS(x)) ) THEN
            DrFP = .TRUE.
        ELSE
            DrFP = .FALSE.
        ENDIF

        !    SOURCE listing of EQ_FP8.f90:
        !       u = DMIN1(DABS(x), DABS(y))
        !       v = DMAX1(DABS(x), DABS(y))
        !       Difx  = NEAREST( u, v )
        !       Dify  = SPACING( v )
        !       Difxy = NEAREST( v-u, v )
        !       IF (DABS(Difxy) < TINY(u)) Difxy = TINY(u)
        !       IF ( x == y ) THEN
        !           DrFP = .TRUE.
        !       ELSEIF ( DABS(Difxy) < 2.D0*DABS(Difx + Dify - DABS(x)) ) THEN
        !           DrFP = .TRUE.
        !       ELSE
        !           DrFP = .FALSE.
        !       ENDIF
Frank

  • jeff
  • Administrator
  • Offline

Re: No Intrinsic Procedure for real equals real?

Frank,

This post is absolutely appropriate for this forum.  Please feel free to post things like this.  I'll try to answer your questions as best I can.

(1) Is this the best way to test for equality?

It really depends on the situation.  If you have a simple calculation that you think the floating point result should be exactly a number, then this routine will work fine.  The problem arises when you're testing a drastically more complex calculation that you think should end up "equal" to a given value.  The errors inherent in floating-point mathematics on a computer can build as calculations use the results of previous calculations.  While this build-up of errors might result in a value that is different from the expected answer by more than machine precision (i.e. the result of Fortran's NEAREST intrinsic), it is still effectively correct.

In my work, I would usually make a subjective judgment based on the problem being solved.  For example, for some temperatures resulting from slow combustion, it was more than sufficient to test for equality within 0.01K.  There was no need to get closer considering the flame was burning near 500-600K.

I normally will just use something like

IF(ABS(x-y) < 0.01) THEN ! Equal

(2) Is there a better approach that works well?

Again, it would depend on what you're trying to do.  You can also think about how complicated you want the check to be.  For example, your equality function is great, but if you're calling it in a loop that iterates over 400,000 points in two arrays every time step, maybe you'll start spending too much time in that calculation.  It might be better to use something more like what I suggested above that only involves subtraction, a single bit modification, and a comparison.

(3) Why does double precision x+factor == y1 (.true.) in my example?

Your value of y1 is the next largest representable value from x.  In your test routine, though, you subtract the two numbers and then find the next largest value (Difxy) from that difference that is representable. 

This Difxy value is then tested against:

2.D0*DABS(Difx + Dify - DABS(x))

This is adding the smaller (x) to the spacing of numbers calculated from the larger (SPACING(y1)) and subtracting the absolute value of x explicitly.  The answer from this calculation should just be SPACING(y1), which should effectively be the same as SPACING(x).  That value is multiplied by two, which I assume is almost like a buffer for further floating-point inaccuracy.  But because Difxy is effectively going to be SPACING(x) anyway, it is indeed smaller than double the same value.

The same test fails for single precision, of course, because Difxy is computed in double precision, whereas the difference is actually going to be much larger because x and y were created in single precision.

I would say that you might want to create a module that implements this same routine in single-, double-, and quad-precision and uses an interface to call the appropriate routine based on precision.  Something like:

INTERFACE EQ_FP
    MODULE PROCEDURE EQ_FP4, EQ_FP8, EQ_FP16
END INTERFACE EQ_FP

Then code could just call EQ_FP, and the interface block would route it to the proper precision test.

Jeff Armstrong
Approximatrix, LLC

Re: No Intrinsic Procedure for real equals real?

Hi Jeff,

Thank you for your good suggestions on improving the equality function for real numbers.

Your suggestions were,
"The answer from this calculation should just be SPACING(y1), which should effectively be the same as SPACING(x)."
and
"you might want to create a module that implements this same routine in single-, double-, and quad-precision and uses an interface to call the appropriate routine based on precision."

Per your suggestions, please find a new and improved routine for testing the equality of two real numbers, based on SPACING and an interface module for the appropriate precision routine, listed below.

KIND regards,
Frank

!--- The EQUALITY MODULE is listed here:

        MODULE EQUAL_FP
        !***************************************************************
        !---    EQ_FP version 1.1
        !---
        !---    Developed by F.W. Perrella, Ph.D.  November 9, 2020
        !---    EQ_FP™ Copyright © 2020 F.W. PERRELLA, Ph.D.
        !---            Silver Spring, Maryland U.S.A.
        !---                All Rights Reserved.
        !***************************************************************
        ! NAME:
        !       EQ_FP
        !
        ! PURPOSE:
        !       Program to determine whether two real numbers are equal.
        !
        ! MODULE:
        !       Type_Kinds: Definition for kind variable types.
        !
        ! ROUTINE:
        !       EQ_FP: Routine to perform equality comparison of
        !              single, double, quadruple precision real numbers.
        !---------------------------------------------------------------
        !
        !    CALLING SEQUENCE:
        !       Result = EQ_FP( x,y )
        !
        !    INPUT ARGUMENTS:
        !       x, y:       Two real numbers to compare.
        !       TYPE:       REAL( KIND = Single, Double, Quadruple )
        !
        !    FUNCTION RESULT:
        !       TYPE:    Logical
        !
        !       RESULT: A returned logical value indicating whether the
        !               x and y inputs are equal within the set precision.
        !               Set to .TRUE. if the real numbers are equal
        !               within the specified tolerance.
        !               Set to .FALSE. if the real numbers are not equal.
        !---------------------------------------------------------------
        USE Type_Kinds

        IMPLICIT NONE

        !--------------------
        !--- Visible function
        !--------------------
        PUBLIC :: EQ_FP

        !----------------------------
        !--- Floating Point Precision
        !----------------------------
        !    REAL( KIND )
        !============================
        !    Single (SP, 4  bytes)
        !    Double (DP, 8  bytes)
        !    Quad   (QP, 16 bytes)
        !----------------------------

        !---------------------
        !--- Procedures
        !---------------------
        INTERFACE EQ_FP
            MODULE PROCEDURE EQ_FP4     !--- Kind(SP)
            MODULE PROCEDURE EQ_FP8     !--- Kind(DP)
            MODULE PROCEDURE EQ_FP16    !--- Kind(QP)
        END INTERFACE EQ_FP

        !---------------------------------------------------------------
        CONTAINS
        !---------------------------------------------------------------

        !---------------------------------------------------------------
        ELEMENTAL FUNCTION EQ_FP4( x,y ) Result( IsEQ )
        !--- USE Type_Kinds
        IMPLICIT NONE

        !--- Define Kind as single precision
        INTEGER, PARAMETER             :: KIND_FP = SP

        LOGICAL                        :: IsEQ
        REAL( KIND_FP ), INTENT( IN )  :: x
        REAL( KIND_FP ), INTENT( IN )  :: y

        !--- Local Variables
        REAL( KIND_FP )                :: DiffR
        REAL( KIND_FP )                :: Dif
        REAL( KIND_FP ), PARAMETER     :: Rel = 2.0_KIND_FP
        INTRINSIC                      :: NEAREST, SPACING, ABS

        Dif   = NEAREST( ABS(y-x),x )
        DiffR = SPACING( ABS(x) )

        !--- Test for equality
        IF ( x == y ) THEN
            IsEQ = .TRUE.
        ELSEIF ( Dif < Rel*DiffR ) THEN
            IsEQ = .TRUE.
        ELSE
            IsEQ = .FALSE.
        ENDIF

        END FUNCTION EQ_FP4

        !---------------------------------------------------------------
        ELEMENTAL FUNCTION EQ_FP8( x,y ) Result( IsEQ )
        !--- USE Type_Kinds
        IMPLICIT NONE

        !--- Define Kind as double precision
        INTEGER, PARAMETER             :: KIND_FP = DP

        LOGICAL                        :: IsEQ
        REAL( KIND_FP ), INTENT( IN )  :: x
        REAL( KIND_FP ), INTENT( IN )  :: y

        !--- Local Variables
        REAL( KIND_FP )                :: DiffR
        REAL( KIND_FP )                :: Dif
        REAL( KIND_FP ), PARAMETER     :: Rel = 2.0_KIND_FP
        INTRINSIC                      :: NEAREST, SPACING, ABS

        Dif   = NEAREST( ABS(y-x),x )
        DiffR = SPACING( ABS(x) )

        !--- Test for equality
        IF ( x == y ) THEN
            IsEQ = .TRUE.
        ELSEIF ( ABS(Dif) < Rel*DiffR ) THEN
            IsEQ = .TRUE.
        ELSE
            IsEQ = .FALSE.
        ENDIF

        END FUNCTION EQ_FP8

        !---------------------------------------------------------------
        ELEMENTAL FUNCTION EQ_FP16( x,y ) Result( IsEQ )
        !--- USE Type_Kinds
        IMPLICIT NONE

        !--- Define Kind as quadruple precision
        INTEGER, PARAMETER             :: KIND_FP = QP

        LOGICAL                        :: IsEQ
        REAL( KIND_FP ), INTENT( IN )  :: x
        REAL( KIND_FP ), INTENT( IN )  :: y

        !--- Local Variables
        REAL( KIND_FP )                :: DiffR
        REAL( KIND_FP )                :: Dif
        REAL( KIND_FP ), PARAMETER     :: Rel = 2.0_KIND_FP
        INTRINSIC                      :: NEAREST, SPACING, ABS

        Dif   = NEAREST( ABS(y-x),x )
        DiffR = SPACING( ABS(x) )

        !--- Test for equality
        IF ( x == y ) THEN
            IsEQ = .TRUE.
        ELSEIF ( ABS(Dif) < Rel*DiffR ) THEN
            IsEQ = .TRUE.
        ELSE
            IsEQ = .FALSE.
        ENDIF

        END FUNCTION EQ_FP16
        !---------------------------------------------------------------
        END MODULE EQUAL_FP
        !---------------------------------------------------------------
       
!--- The MAIN PROGRAM is listed here:

        PROGRAM Test_EQUAL_FP
        !---------------------------------------------------------------
        ! NAME:
        !       Test_EQ_FP
        ! PURPOSE:
        !       Program to test the routine EQ_FP(x,y).
        ! MODULES:
        !       Type_Kinds: Definitions for kinds of variable types.
        !       EQ_FP:   Performs equality comparison of real numbers.
        !---------------------------------------------------------------
        USE Type_Kinds
        USE EQUAL_FP

        IMPLICIT NONE
        !-------------------------------
        !    Floating Point Precision
        !-------------------------------
        !    REAL( KIND )
        !===============================
        !    Single (SP, 4  bytes)
        !    Double (DP, 8  bytes)
        !    Quad   (QP, 16 bytes)
        !-------------------------------

        !-----------
        ! Parameters
        !-----------
        CHARACTER( * ), PARAMETER :: Title = 'Test EQ_FP(x,y):'
        CHARACTER( * ), PARAMETER :: cLine = REPEAT("-",50)
        CHARACTER( * ), PARAMETER :: cDot = REPEAT(".",50)

        !---------------------------------
        ! Numbers to test for precision
        !---------------------------------
        INTEGER, PARAMETER :: NUMBERS = 5
        REAL( SP ), PARAMETER, DIMENSION( NUMBERS ) :: SINGLE_NUMBER = &
        (/ 1.000000000000000e-32_SP,    &
        1.234567890123456e-16_SP,       &
        1.000000000000000e+00_SP,       &
        1.234567890123456e+16_SP,       &
        1.234567890123456e+32_SP /)

        REAL( DP ), PARAMETER, DIMENSION( NUMBERS ) :: DOUBLE_NUMBER = &
        (/ 1.000000000000000e-32_DP,    &
        1.234567890123456e-16_DP,       &
        1.000000000000000e+00_DP,       &
        1.234567890123456e+16_DP,       &
        1.234567890123456e+32_DP /)

        REAL( QP ), PARAMETER, DIMENSION( NUMBERS ) :: QUAD_NUMBER = &
        (/ 1.000000000000000e-32_QP,  &
        1.234567890123456e-16_QP,     &
        1.000000000000000e+00_QP,     &
        1.234567890123456e+16_QP,     &
        1.234567890123456e+32_QP /)

        !----------
        ! Variables
        !----------
        CHARACTER(len=256) :: cend
        INTEGER :: i
        REAL( SP ) ::  x,  y1,  y2,  y3,  y4,  x0,  y0, es
        REAL( DP ) :: xd, yd1, yd2, yd3, yd4, xd0, yd0, ed
        REAL( QP ) :: xq, yq1, yq2, yq3, yq4, xq0, yq0, eq

        !-------------------------------
        ! Test equality of x and y reals
        !-------------------------------
        WRITE( *, '(/5x,a)') Title
        WRITE( *, '(5x,a)')  cLine
        DO i = 1, NUMBERS
            WRITE( *, '(5x, "NUMBER ", i2, " of ", i2 )' ) i, NUMBERS
            !----------------------
            ! Single precision test
            !----------------------
            x = SINGLE_NUMBER(i)
            y1 = NEAREST( x, 1.0_SP )
            y2 = y1 - SPACING( x )
            y3 = NEAREST( x, -1.0_SP )
            y4 = y3 + SPACING( x )

            !----------------
            ! Error variables
            !----------------
            es = x * 5.e-08_SP
            x = x + es

            !--- Test single precision function
            WRITE( *, '( /5x, "SINGLE PRECISION: x  = ", es20.13, &
            &/5x, "y1 = ", es20.13, 2x, ":  NEAREST( x, 1.0 )", &
            &/5x, "y2 = ", es20.13, 2x, ":  y1 - SPACING( x )", &
            &/5x, "y3 = ", es20.13, 2x, ":  NEAREST( x,-1.0 )", &
            &/5x, "y4 = ", es20.13, 2x, ":  y3 + SPACING( x )" )' ) &
            x, y1, y2, y3, y4

            x0 = x; y0 = y1
            WRITE( *, '(/5x,a,2x,L6)') '(x==y1): X=Y?', (x0==y0)
            x0 = x; y0 = y2
            WRITE( *, '(5x,a,2x,L6)') '(x==y1): X=Y?', (x0==y0)
            x0 = x; y0 = y3
            WRITE( *, '(5x,a,2x,L6)') '(x==y1): X=Y?', (x0==y0)
            x0 = x; y0 = y4
            WRITE( *, '(5x,a,2x,L6)') '(x==y1): X=Y?', (x0==y0)

            x0 = x; y0 = y1
            WRITE( *, '(/5x,a,2x,L6)') 'EQ_FP4(x,y1): X=Y?', EQ_FP(x0,y0)
            x0 = x; y0 = y2
            WRITE( *, '( 5x,a,2x,L6)') 'EQ_FP4(x,y2): X=Y?', EQ_FP(x0,y0)
            x0 = x; y0 = y3
            WRITE( *, '( 5x,a,2x,L6)') 'EQ_FP4(x,y3): X=Y?', EQ_FP(x0,y0)
            x0 = x; y0 = y4
            WRITE( *, '( 5x,a,2x,L6)') 'EQ_FP4(x,y4): X=Y?', EQ_FP(x0,y0)

            !----------------------
            ! Double precision test
            !----------------------
            xd = DOUBLE_NUMBER(i)
            yd1 = NEAREST( xd, 1.0_DP )
            yd2 = yd1 - SPACING( xd )
            yd3 = NEAREST( xd, -1.0_DP )
            yd4 = yd3 + SPACING( xd )

            !----------------
            ! Error variables
            !----------------
            ed = xd * 5.e-16_DP
            xd = xd + ed

            !--- Test stand alone double precision function routines
            WRITE( *, '(/5x,a)')  cDot
            WRITE( *, '( /5x, "DOUBLE PRECISION: x  = ", es20.13, &
            &/5x, "y1 = ", es20.13, 2x, ":  NEAREST( x, 1.0 )", &
            &/5x, "y2 = ", es20.13, 2x, ":  y1 - SPACING( x )", &
            &/5x, "y3 = ", es20.13, 2x, ":  NEAREST( x,-1.0 )", &
            &/5x, "y4 = ", es20.13, 2x, ":  y3 + SPACING( x )" )' ) &
            xd, yd1, yd2, yd3, yd4

            xd0 = xd; yd0 = yd1
            WRITE( *, '(/5x,a,2x,L6)') '(x==y1): X=Y?', (x0==y0)
            xd0 = xd; yd0 = yd2
            WRITE( *, '(5x,a,2x,L6)') '(x==y1): X=Y?', (x0==y0)
            xd0 = xd; yd0 = yd3
            WRITE( *, '(5x,a,2x,L6)') '(x==y1): X=Y?', (x0==y0)
            xd0 = xd; yd0 = yd4
            WRITE( *, '(5x,a,2x,L6)') '(x==y1): X=Y?', (x0==y0)

            xd0 = xd; yd0 = yd1
            WRITE( *, '(/5x,a,2x,L6)') 'EQ_FP8(x,y1): X=Y?', EQ_FP(x0,y0)
            xd0 = xd; yd0 = yd2
            WRITE( *, '(5x,a,2x,L6)') 'EQ_FP8(x,y2): X=Y?', EQ_FP(x0,y0)
            xd0 = xd; yd0 = yd3
            WRITE( *, '(5x,a,2x,L6)') 'EQ_FP8(x,y3): X=Y?', EQ_FP(x0,y0)
            xd0 = xd; yd0 = yd4
            WRITE( *, '(5x,a,2x,L6)') 'EQ_FP8(x,y4): X=Y?', EQ_FP(x0,y0)

            !-------------------------
            ! Quadruple precision test
            !-------------------------
            xq = QUAD_NUMBER(i)
            yq1 = NEAREST( xq, 1.0_QP )
            yq2 = yq1 - SPACING( xq )
            yq3 = NEAREST( xq, -1.0_QP )
            yq4 = yq3 + SPACING( xq )

            !----------------
            ! Error variables
            !----------------
            eq = xq * 5.e-32_QP
            xq = xq + eq

            !--- Test quadruple precision function
            WRITE( *, '(/5x,a)')  cDot
            WRITE( *, '( /5x, "QUADRULE PRECISION: x  = ", es20.13, &
            &/5x, "y1 = ", es20.13, 2x, ":  NEAREST( x, 1.0 )", &
            &/5x, "y2 = ", es20.13, 2x, ":  y1 - SPACING( x )", &
            &/5x, "y3 = ", es20.13, 2x, ":  NEAREST( x,-1.0 )", &
            &/5x, "y4 = ", es20.13, 2x, ":  y3 + SPACING( x )" )' ) &
            xq, yq1, yq2, yq3, yq4

            xq0 = xq; yq0 = yq1
            WRITE( *, '(/5x,a,2x,L6)') '(x==y1): X=Y?', (x0==y0)
            xq0 = xq; yq0 = yq2
            WRITE( *, '(5x,a,2x,L6)') '(x==y1): X=Y?', (x0==y0)
            xq0 = xq; yq0 = yq3
            WRITE( *, '(5x,a,2x,L6)') '(x==y1): X=Y?', (x0==y0)
            xq0 = xq; yq0 = yq4
            WRITE( *, '(5x,a,2x,L6)') '(x==y1): X=Y?', (x0==y0)

            xq0 = xq; yq0 = yq1
            WRITE( *, '(/5x,a,2x,L6)') 'EQ_FP16(x,y1): X=Y?', EQ_FP(x0,y0)
            xq0 = xq; yq0 = yq2
            WRITE( *, '(5x,a,2x,L6)') 'EQ_FP16(x,y2): X=Y?', EQ_FP(x0,y0)
            xq0 = xq; yq0 = yq3
            WRITE( *, '(5x,a,2x,L6)') 'EQ_FP16(x,y3): X=Y?', EQ_FP(x0,y0)
            xq0 = xq; yq0 = yq4
            WRITE( *, '(5x,a,2x,L6)') 'EQ_FP16(x,y4): X=Y?', EQ_FP(x0,y0)

            WRITE( *, '(/5x,a)')  cLine
        END DO

        WRITE( *, '( /10x, "Press <ENTER> to continue..." )' )
        READ( *, '(a)' ) cend

        END PROGRAM Test_EQUAL_FP
        !---------------------------------------------------------------

!--- The KIND TYPE MODULE is listed here:

    MODULE Type_Kinds
        IMPLICIT NONE
        PRIVATE
        !-------------------------
        ! Floating Point Precision
        !-------------------------
        !    REAL( KIND )
        !=========================
        !    Single (SP, 4  bytes)
        !    Double (DP, 8  bytes)
        !    Quad   (QP, 16 bytes)
        !-------------------------
        !--- Single precision
        INTEGER, PARAMETER :: &
        Single = SELECTED_REAL_KIND(6)

        !--- Double precision
        INTEGER, PARAMETER :: &
        Double = SELECTED_REAL_KIND(2*precision(1.0_Single))

        !--- Quadruple precision
        INTEGER, PARAMETER :: &
        Quadr = SELECTED_REAL_KIND(2*precision(1.0_Double))
        INTEGER, PARAMETER :: &
        Quad = (((1 + SIGN(1,Quadr)) / 2) * Quadr) &
        + (((1 - SIGN(1,Quadr)) / 2) * Double)

        INTEGER, PARAMETER, PUBLIC:: SP = Single, DP = Double, QP = Quad

    END MODULE Type_Kinds
   
!--- The END of the program EQ_FP(x,y)
! -----------------------------------------------------------------------------------------------------