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)
! -----------------------------------------------------------------------------------------------------
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_FPThen code could just call EQ_FP, and the interface block would route it to the proper precision test.
]]>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
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
!---------------------------------------------------------------
!---------------------------------------------------------------