! This file: http://ftp.aset.psu.edu/pub/ger/fortran/hdk/SpacingDemo.f90
!
  Program SpacingDemo
!  Purpose: To demonstrate use of Spacing and Epsilon.
!
!  by Paul Van Delst, James Van Buskirk, and James Giles.
!   Posted to Comp.Lang.Fortran on 16 Sep 2005.
!
!  Source: http://ftp.cac.psu.edu/pub/ger/fortran/hdk/SpacingDemo.f90
!  Also see: http://ftp.cac.psu.edu/pub/ger/fortran/hdk/Eps.f90
!
! Compare two real numbers x and y.
! ABS( x - y ) < ( ULP * SPACING( MAX(ABS(x),ABS(y)) ) )
!
!       The test performed is
!
!         ABS( x - y ) < ( ULP * SPACING( MAX(ABS(x),ABS(y)) ) )
!
!       If the result is .TRUE., the numbers are considered equal.
!
!       The intrinsic function SPACING(x) returns the absolute spacing
!       of numbers near the value of x,
!
!                      {     EXPONENT(x)-DIGITS(x)
!                      {  2.0                        for x /= 0
!         SPACING(x) = {
!                      {
!                      {  TINY(x)                    for x == 0
!
!       The ULP optional argument scales the comparison.
!
!       James Van Buskirk and James Giles suggested this method for
!       floating-point comparisons in the comp.lang.fortran newsgroup.
!
! Where "ULP" is a user supplied scale factor (for those case where your
! numbers are close, but not that close). Definition of ULP:
!
! OPTIONAL INPUT ARGUMENTS:
!  ULP: Unit of data precision. The acronym stands for "unit in
!       the last place," the smallest possible increment or decrement
!       that can be made using a machine"s floating point arithmetic.
!       A 0.5 ulp maximum error is the best you could hope for, since
!       this corresponds to always rounding to the nearest representable
!       floating-point number. Value must be positive - if a negative
!       value is supplied, the absolute value is used.
!       If not specified, the default value is 1.

   implicit none
   double precision  x, a, b, fuzz, ulp
   ulp = 10.0D0;  fuzz=EPSILON(x)*ulp
   print *, " "
   print *, "ulp=", ulp,"  fuzz=",fuzz

   if (Spacing(1.0D0) == Epsilon(x)) then
     print *, "Spacing(1.0D0) = Epsilon = ",Spacing(1.0D0)
   else
     print *, "Spacing(1.0D0) /= Epsilon"
   endif
   if (Spacing(0.0D0) == Tiny(x)) then
     print *, "Spacing(0.0D0) = Tiny=",Tiny(x)
   else
     print *, "Spacing(0.0D0) /= Tiny"
   endif
   print *, " "

!  Define nearly equal quantities, a, b
   a=0.1234567898765443D+20
   b=0.1234567898765441D+20
   print *, " "
   print *, "a,b=",a,b
   print *, "a-b=",a-b
   print *, "Spacing(a)=",Spacing(a)
   print *, " "

! Do absolute comparison
   if (a == b) then
     print *, "Absolute comparison a==b"
   else
     print *, "Absolute comparison a/=b"
   endif
!  Do comparisons using Spacing and Epsilon intrinsic procedures.
   if(abs(a-b) <= ulp*Spacing( Max(Abs(a),Abs(b)))) then
     print *, "Spacing comparison: Equal."
   else
     print *, "Spacing comparison: Not Equal."
   endif

   if (abs(a-b) <= MAX(abs(a),abs(b))*fuzz) then
     print *, "Epsilon comparison: Equal."
   else
     print *, "Epsilon comparison: Not Equal."
   endif

!  Output:
!
!      ulp= 10.00000000000000   fuzz= 2.220446049250313E-15
!      Spacing(1.0D0) = Epsilon =  2.220446049250313E-16
!      Spacing(0.0D0) = Tiny= 2.225073858507201E-308
!
!
!      a,b= 1.234567898765443E+19 1.234567898765441E+19
!      a-b= 20480.00000000000
!      Spacing(a)= 2048.000000000000
!
!      Absolute comparison a/=b
!      Spacing comparison: Equal.
!      Epsilon comparison: Equal.

 end program