Compare_Float_Numbers.f90

Go to the documentation of this file.

!
! Compare_Float_Numbers
!
! Module containing routines to perform equality and relational
! comparisons on floating point numbers.
!
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 01-Apr-2003
!                       paul.vandelst@noaa.gov
!

MODULE compare_float_numbers


  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module usage
  USE type_kinds, ONLY: single, double
  ! Disable all implicit typing
  IMPLICIT NONE


  ! ------------
  ! Visibilities
  ! ------------
  PRIVATE
  ! Parameters
  PUBLIC :: default_n_sigfig
  ! Operators
  PUBLIC :: OPERATOR (.equalto.)
  PUBLIC :: OPERATOR (.greaterthan.)
  PUBLIC :: OPERATOR (.lessthan.)
  ! Procedures
  PUBLIC :: compare_float
  PUBLIC :: tolerance
  PUBLIC :: compares_within_tolerance


  ! ---------------------
  ! Procedure overloading
  ! ---------------------

  INTERFACE compare_float
    MODULE PROCEDURE compare_real_single
    MODULE PROCEDURE compare_real_double
    MODULE PROCEDURE compare_complex_single
    MODULE PROCEDURE compare_complex_double
  END INTERFACE compare_float

  INTERFACE OPERATOR (.equalto.)
    MODULE PROCEDURE equalto_real_single
    MODULE PROCEDURE equalto_real_double
    MODULE PROCEDURE equalto_complex_single
    MODULE PROCEDURE equalto_complex_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.)

  INTERFACE tolerance
    MODULE PROCEDURE tolerance_real_single
    MODULE PROCEDURE tolerance_real_double
    MODULE PROCEDURE tolerance_complex_single
    MODULE PROCEDURE tolerance_complex_double
  END INTERFACE tolerance

  INTERFACE compares_within_tolerance
    MODULE PROCEDURE cwt_real_single
    MODULE PROCEDURE cwt_real_double
    MODULE PROCEDURE cwt_complex_single
    MODULE PROCEDURE cwt_complex_double
  END INTERFACE compares_within_tolerance

  ! -----------------
  ! Module parameters
  ! -----------------
  ! Module Version Id
  CHARACTER(*), PARAMETER :: module_version_id = &
    '$Id: Compare_Float_Numbers.f90 60152 2015-08-13 19:19:13Z paul.vandelst@noaa.gov $'
  ! Numeric literals
  REAL(Single), PARAMETER :: sp_zero = 0.0_single
  REAL(Double), PARAMETER :: dp_zero = 0.0_double
  REAL(Single), PARAMETER :: sp_one = 1.0_single
  REAL(Double), PARAMETER :: dp_one = 1.0_double
  REAL(Single), PARAMETER :: sp_ten = 10.0_single
  REAL(Double), PARAMETER :: dp_ten = 10.0_double
  REAL(Single), PARAMETER :: sp_hundred = 100.0_single
  REAL(Double), PARAMETER :: dp_hundred = 100.0_double
  REAL(Single), PARAMETER :: sp_compare_cutoff = 1.0e-15_single
  REAL(Double), PARAMETER :: dp_compare_cutoff = 1.0e-15_double
  ! Default number of significant figures
  INTEGER, PARAMETER :: default_n_sigfig = 6


CONTAINS


!----------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!       .EqualTo.
!
! PURPOSE:
!       Relational operator to test the equality of REAL operands.
!
! CALLING SEQUENCE:
!       IF ( x .EqualTo. y ) THEN
!         .....
!       END IF
!
! OPERANDS:
!       x, y:        Two congruent floating point data objects to compare.
!                    UNITS:      N/A
!                    TYPE:       REAL(Single),    REAL(Double)
!                                COMPLEX(Single), COMPLEX(Double)
!                    DIMENSION:  Scalar, or any allowed rank array.
!
! OPERATOR RESULT:
!       (x .EqualTo. y)    The result is a logical value indicating whether
!                          the operands are equal to within numerical precision
!                          UNITS:      N/A
!                          TYPE:       LOGICAL
!                          DIMENSION:  Same as operands.
!
! PROCEDURE:
!       The test performed is
!
!         ABS( x - y ) < SPACING( MAX(ABS(x),ABS(y)) )
!
!       If the result is .TRUE., the numbers are considered equal. For complex
!       input the test is applied separately to the real and imaginary parts.
!:sdoc-:
!----------------------------------------------------------------------------------

  ELEMENTAL FUNCTION equalto_real_single( x, y ) RESULT( EqualTo )
    REAL(Single), INTENT(IN)  :: x, y
    LOGICAL :: equalto
    equalto = abs(x-y) < spacing( max(abs(x),abs(y)) )
  END FUNCTION equalto_real_single

  ELEMENTAL FUNCTION equalto_real_double( x, y ) RESULT( EqualTo )
    REAL(Double), INTENT(IN)  :: x, y
    LOGICAL :: equalto
    equalto = abs(x-y) < spacing( max(abs(x),abs(y)) )
  END FUNCTION equalto_real_double

  ELEMENTAL FUNCTION equalto_complex_single( x, y ) RESULT( EqualTo )
    COMPLEX(Single), INTENT(IN)  :: x, y
    LOGICAL :: equalto
    REAL(Single) :: rx, ix
    REAL(Single) :: ry, iy
    rx = REAL(x,single); ix = aimag(x)
    ry = REAL(y,single); iy = aimag(y)
    equalto = equalto_real_single( rx, ry ) .AND. equalto_real_single( ix, iy )
  END FUNCTION equalto_complex_single

  ELEMENTAL FUNCTION equalto_complex_double( x, y ) RESULT( EqualTo )
    COMPLEX(Double), INTENT(IN)  :: x, y
    LOGICAL :: equalto
    REAL(Double) :: rx, ix
    REAL(Double) :: ry, iy
    rx = REAL(x,double); ix = aimag(x)
    ry = REAL(y,double); iy = aimag(y)
    equalto = equalto_real_double( rx, ry ) .AND. equalto_real_double( ix, iy )
  END FUNCTION equalto_complex_double


!----------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!       .GreaterThan.
!
! PURPOSE:
!       Relational operator to test if one REAL operand is greater than another.
!
! CALLING SEQUENCE:
!       IF ( x .GreaterThan. y ) THEN
!         .....
!       END IF
!
! OPERANDS:
!       x, y:        Two congruent floating point data objects to compare.
!                    UNITS:      N/A
!                    TYPE:       REAL(Single)   [ == default real]
!                                  OR
!                                REAL(Double)
!                    DIMENSION:  Scalar, or any allowed rank array.
!
! OPERATOR RESULT:
!       (x .GreaterThan. y)    The result is a logical value indicating whether
!                              the operand x is greater than y by more than
!                              the spacing between representable floating point
!                              numbers.
!                              UNITS:      N/A
!                              TYPE:       LOGICAL
!                              DIMENSION:  Same as operands.
!
! PROCEDURE:
!       The test performed is
!
!         ( x - y ) >= SPACING( MAX(ABS(x),ABS(y)) )
!
!       If the result is .TRUE., x is considered greater than y.
!:sdoc-:
!----------------------------------------------------------------------------------

  ELEMENTAL FUNCTION is_greater_than_single( x, y ) RESULT ( Greater_Than )
    REAL(Single), 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_single


  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


!----------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!       .LessThan.
!
! PURPOSE:
!       Relational operator to test if one REAL operand is less than another.
!
! CALLING SEQUENCE:
!       IF ( x .LessThan. y ) THEN
!         .....
!       END IF
!
! OPERANDS:
!       x, y:        Two congruent floating point data objects to compare.
!                    UNITS:      N/A
!                    TYPE:       REAL(Single)   [ == default real]
!                                  OR
!                                REAL(Double)
!                    DIMENSION:  Scalar, or any allowed rank array.
!
! OPERATOR RESULT:
!       (x .LessThan. y)    The result is a logical value indicating whether
!                           the operand x is less than y by more than the
!                           spacing between representable floating point
!                           numbers.
!                           UNITS:      N/A
!                           TYPE:       LOGICAL
!                           DIMENSION:  Same as operands.
!
! PROCEDURE:
!       The test performed is
!
!         ( y - x ) >= SPACING( MAX(ABS(x),ABS(y)) )
!
!       If the result is .TRUE., x is considered less than y.
!:sdoc-:
!----------------------------------------------------------------------------------

  ELEMENTAL FUNCTION is_less_than_single( x, y ) RESULT ( Less_Than )
    REAL(Single), 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_single


  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


!----------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!       Compare_Float
!
! PURPOSE:
!       Function to compare floating point scalars and arrays with adjustible
!       precision tolerance.
!
! CALLING SEQUENCE:
!       Result = Compare_Float( x, y,            &  ! Input
!                               ULP    =ULP    , &  ! Optional input
!                               Percent=Percent  )  ! Optional input
!
! INPUT ARGUMENTS:
!       x, y:        Two congruent floating point data objects to compare.
!                    UNITS:      N/A
!                    TYPE:       REAL(Single)   [ == default real]
!                                  OR
!                                REAL(Double)
!                                  OR
!                                COMPLEX(Single)
!                                  OR
!                                COMPLEX(Double)
!                    DIMENSION:  Scalar, or any allowed rank array.
!                    ATTRIBUTES: INTENT(IN)
!
! 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.
!                    This argument is ignored if the Percent optioanl argument is specifed.
!                    UNITS:      N/A
!                    TYPE:       INTEGER
!                    DIMENSION:  Scalar
!                    ATTRIBUTES: OPTIONAL, INTENT(IN)
!
!       Percent:     Specify a percentage difference value to use in comparing
!                    the numbers rather than testing within some numerical
!                    limit. The ULP argument is ignored if this argument is
!                    specified.
!                    UNITS:      N/A
!                    TYPE:       REAL(Single)  for REAL(Single) or COMPLEX(Single) x,y
!                                  OR
!                                REAL(Double)  for REAL(Double) or COMPLEX(Double) x,y
!                    DIMENSION:  Scalar
!                    ATTRIBUTES: OPTIONAL, INTENT(IN)
!
! FUNCTION RESULT:
!       Result:      The return value is a logical value indicating whether
!                    the inputs are equal (to within the required precision)
!                    .TRUE.  - if the floating point numbers are equal to
!                              within the specified tolerance. 
!                    .FALSE. - if the floating point numbers are different.
!                    UNITS:      N/A
!                    TYPE:       LOGICAL
!                    DIMENSION:  Scalar
!
! PROCEDURE:
!       ULP Test
!       --------
!       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.
!
!
!       Percent Test
!       ------------
!       The test performed is
!
!         100.0 * ABS((x-y)/x) < Percent
!
!       If the result is .TRUE., the numbers are considered equal.
!
!
!       For complex numbers, the same test is applied to both the real and
!       imaginary parts and each result is ANDed.
!:sdoc-:
!----------------------------------------------------------------------------------

  ELEMENTAL FUNCTION compare_real_single( x, y, ULP, Percent ) RESULT( Compare )
    ! Arguments
    REAL(Single),           INTENT(IN) :: x
    REAL(Single),           INTENT(IN) :: y
    INTEGER     , OPTIONAL, INTENT(IN) :: ulp
    REAL(Single), OPTIONAL, INTENT(IN) :: percent
    ! Function result
    LOGICAL :: compare
    ! Local variables
    LOGICAL      :: ulp_test
    REAL(Single) :: rel
    
    ! Set up
    ! ------
    ulp_test = .true.
    IF ( PRESENT(ulp) ) THEN
      rel = REAL(ABS(ULP), single)
    ELSE
      rel = sp_one
    END IF
    IF ( PRESENT(percent) ) THEN
      ulp_test = .false.
      ! Test for zero x (elementals can't be recursive)
      IF ( abs(x) < ( spacing( max(abs(x),sp_zero) ) ) ) ulp_test = .true.
    END IF
    
    ! Compare the numbers
    ! -------------------
    IF ( ulp_test ) THEN
      compare = abs(x-y) < ( rel * spacing( max(abs(x),abs(y)) ) )
    ELSE
      compare = sp_hundred*abs((x-y)/x) < percent
    END IF
  END FUNCTION compare_real_single


  ELEMENTAL FUNCTION compare_real_double( x, y, ULP, Percent ) RESULT( Compare )
    ! Arguments
    REAL(Double),           INTENT(IN) :: x
    REAL(Double),           INTENT(IN) :: y
    INTEGER     , OPTIONAL, INTENT(IN) :: ulp
    REAL(Double), OPTIONAL, INTENT(IN) :: percent
    ! Function result
    LOGICAL :: compare
    ! Local variables
    LOGICAL      :: ulp_test
    REAL(Double) :: rel
    
    ! Set up
    ! ------
    ulp_test = .true.
    IF ( PRESENT(ulp) ) THEN
      rel = REAL(ABS(ULP), double)
    ELSE
      rel = dp_one
    END IF
    IF ( PRESENT(percent) ) THEN
      ulp_test = .false.
      ! Test for zero x (elementals can't be recursive)
      IF ( abs(x) < ( spacing( max(abs(x),dp_zero) ) ) ) ulp_test = .true.
    END IF
    
    ! Compare the numbers
    ! -------------------
    IF ( ulp_test ) THEN
      compare = abs(x-y) < ( rel * spacing( max(abs(x),abs(y)) ) )
    ELSE
      compare = dp_hundred*abs((x-y)/x) < percent
    END IF
  END FUNCTION compare_real_double


  ELEMENTAL FUNCTION compare_complex_single( x, y, ULP, Percent ) RESULT( Compare )
    ! Arguments
    COMPLEX(Single),           INTENT(IN) :: x
    COMPLEX(Single),           INTENT(IN) :: y
    INTEGER        , OPTIONAL, INTENT(IN) :: ulp
    REAL(Single)   , OPTIONAL, INTENT(IN) :: percent
    ! Function result
    LOGICAL :: compare
    ! Local variables
    REAL(Single) :: xr, xi
    REAL(Single) :: yr, yi
    
    ! Separate real and complex parts
    ! -------------------------------
    xr=REAL(x,single); xi=aimag(x)
    yr=REAL(y,single); yi=aimag(y)
    
    ! Compare each part separately
    ! ----------------------------
    compare = compare_real_single(xr,yr,ulp=ulp,percent=percent) .AND. &
              compare_real_single(xi,yi,ulp=ulp,percent=percent)
  END FUNCTION compare_complex_single


  ELEMENTAL FUNCTION compare_complex_double( x, y, ULP, Percent ) RESULT( Compare )
    ! Arguments
    COMPLEX(Double),           INTENT(IN) :: x
    COMPLEX(Double),           INTENT(IN) :: y
    INTEGER        , OPTIONAL, INTENT(IN) :: ulp
    REAL(Double)   , OPTIONAL, INTENT(IN) :: percent
    ! Function result
    LOGICAL :: compare
    ! Local variables
    REAL(Double) :: xr, xi
    REAL(Double) :: yr, yi
    
    ! Separate real and complex parts
    ! -------------------------------
    xr=REAL(x,double); xi=aimag(x)
    yr=REAL(y,double); yi=aimag(y)
    
    ! Compare each part separately
    ! ----------------------------
    compare = compare_real_double(xr,yr,ulp=ulp,percent=percent) .AND. &
              compare_real_double(xi,yi,ulp=ulp,percent=percent)
  END FUNCTION compare_complex_double


!----------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!       Tolerance
!
! PURPOSE:
!       Elemental function to compute a tolerance value for a given input for a
!       specified number of significant figures.
!
! CALLING SEQUENCE:
!       Result = Tolerance( x, n )
!
! INPUT ARGUMENTS:
!       x:           Floating point value for which a tolerance value is required.
!                    UNITS:      N/A
!                    TYPE:       REAL(Single)   [ == default real]
!                                  OR
!                                REAL(Double)
!                                  OR
!                                COMPLEX(Single)
!                                  OR
!                                COMPLEX(Double)
!                    DIMENSION:  Scalar or any rank array.
!                    ATTRIBUTES: INTENT(IN)
!
!       n:           The approximate number of significant figures for which the 
!                    tolerance is required.
!                    UNITS:      N/A
!                    TYPE:       INTEGER
!                    DIMENSION:  Scalar or same as input x.
!                    ATTRIBUTES: INTENT(IN)
!
! FUNCTION RESULT:
!       Result:      The return value is a tolerance value that can be used to
!                    compare two numbers.
!                    UNITS:      N/A
!                    TYPE:       Same as input x.
!                    DIMENSION:  Same as input x.
!:sdoc-:
!----------------------------------------------------------------------------------

  ELEMENTAL FUNCTION tolerance_real_single(x,n) RESULT( Tolerance )
    REAL(Single), INTENT(IN) :: x
    INTEGER     , INTENT(IN) :: n
    REAL(Single) :: tolerance
    INTEGER :: e
    IF (abs(x) > sp_zero) THEN
      e = floor(log10(abs(x))) - n
      tolerance = sp_ten**e
    ELSE
      tolerance = sp_one
    END IF
  END FUNCTION tolerance_real_single
  
  ELEMENTAL FUNCTION tolerance_real_double(x,n) RESULT( Tolerance )
    REAL(Double), INTENT(IN) :: x
    INTEGER,      INTENT(IN) :: n
    REAL(Double) :: tolerance
    INTEGER :: e
    IF (abs(x) > dp_zero) THEN
      e = floor(log10(abs(x))) - n
      tolerance = dp_ten**e
    ELSE
      tolerance = dp_one
    END IF
  END FUNCTION tolerance_real_double
  
  ELEMENTAL FUNCTION tolerance_complex_single(x,n) RESULT( Tolerance )
    COMPLEX(Single), INTENT(IN) :: x
    INTEGER,         INTENT(IN) :: n
    COMPLEX(Single) :: tolerance
    REAL(Single) :: tr, ti
    tr = tolerance_real_single(REAL(x,Single),n)
    ti = tolerance_real_single(aimag(x),n)
    tolerance = cmplx(tr,ti,single)
  END FUNCTION tolerance_complex_single
  
  ELEMENTAL FUNCTION tolerance_complex_double(x,n) RESULT( Tolerance )
    COMPLEX(Double), INTENT(IN) :: x
    INTEGER,         INTENT(IN) :: n
    COMPLEX(Double) :: tolerance
    REAL(Double) :: tr, ti
    tr = tolerance_real_double(REAL(x,Double),n)
    ti = tolerance_real_double(aimag(x),n)
    tolerance = cmplx(tr,ti,double)
  END FUNCTION tolerance_complex_double


!----------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!       Compares_Within_Tolerance
!
! PURPOSE:
!       Elemental function to determine if two values are comparable withing
!       a given tolerance determined by the number of significant figures
!       used in the comparison.
!
! CALLING SEQUENCE:
!       Result = Compare_Within_Tolerance( x, y, n, cutoff=cutoff )
!
! INPUTS:
!       x, y:        Floating point values to be compared.
!                    UNITS:      N/A
!                    TYPE:       REAL(Single)   [ == default real]
!                                  OR
!                                REAL(Double)
!                                  OR
!                                COMPLEX(Single)
!                                  OR
!                                COMPLEX(Double)
!                    DIMENSION:  Scalar or any rank array.
!                    ATTRIBUTES: INTENT(IN)
!
!       n:           The approximate number of significant figures for which the 
!                    tolerance is required.
!                    UNITS:      N/A
!                    TYPE:       INTEGER
!                    DIMENSION:  Scalar or same as input x, y.
!                    ATTRIBUTES: INTENT(IN)
!
! OPTIONAL INPUTS:
!       cutoff:      Floating point value below which the comparison is not
!                    performed. In this case, the function result will be .TRUE.
!                    If not specified, the default value is 1.0e-15 for real
!                    input, or (1.0e-15,1.0e-15) for complex input.
!                    UNITS:      N/A
!                    TYPE:       Same as input x
!                    DIMENSION:  Scalar or same as input x, y.
!                    ATTRIBUTES: INTENT(IN), OPTIONAL
!
! FUNCTION RESULT:
!       Result:      The return value is a logical value indicating if the 
!                    comparison was successful or not.
!                    If .TRUE. , the two numbers compare within the prescribed
!                                tolerance, or
!                       .FALSE., they do not.
!                    UNITS:      N/A
!                    TYPE:       LOGICAL
!                    DIMENSION:  Same as input x, y.
!:sdoc-:
!----------------------------------------------------------------------------------

  ELEMENTAL FUNCTION cwt_real_single(x,y,n,cutoff) RESULT(is_comparable)
    REAL(Single),           INTENT(IN) :: x, y
    INTEGER,                INTENT(IN) :: n
    REAL(Single), OPTIONAL, INTENT(IN) :: cutoff
    LOGICAL :: is_comparable
    REAL(Single) :: c
    IF ( PRESENT(cutoff) ) THEN
      c = cutoff
    ELSE
      c = sp_compare_cutoff
    END IF
    is_comparable = .true.
    IF ( abs(x) > c .OR. abs(y) > c ) is_comparable = abs(x-y) < tolerance(x,n)
  END FUNCTION cwt_real_single


  ELEMENTAL FUNCTION cwt_real_double(x,y,n,cutoff) RESULT(is_comparable)
    REAL(Double),           INTENT(IN) :: x, y
    INTEGER,                INTENT(IN) :: n
    REAL(Double), OPTIONAL, INTENT(IN) :: cutoff
    LOGICAL :: is_comparable
    REAL(Double) :: c
    IF ( PRESENT(cutoff) ) THEN
      c = cutoff
    ELSE
      c = dp_compare_cutoff
    END IF
    is_comparable = .true.
    IF ( abs(x) > c .OR. abs(y) > c ) is_comparable = abs(x-y) < tolerance(x,n)
  END FUNCTION cwt_real_double


  ELEMENTAL FUNCTION cwt_complex_single(x,y,n,cutoff) RESULT(is_comparable)
    COMPLEX(Single),           INTENT(IN) :: x, y
    INTEGER,                   INTENT(IN) :: n
    COMPLEX(Single), OPTIONAL, INTENT(IN) :: cutoff
    LOGICAL :: is_comparable
    COMPLEX(Single) :: c
    IF ( PRESENT(cutoff) ) THEN
      c = cutoff
    ELSE
      c = cmplx(sp_compare_cutoff,sp_compare_cutoff,single)
    END IF
    is_comparable = cwt_real_single(REAL(x,Single),REAL(y,Single),n,cutoff=real(c,single) ) .AND. &
                    cwt_real_single(aimag(x),aimag(y),n,cutoff=aimag(c))
  END FUNCTION cwt_complex_single


  ELEMENTAL FUNCTION cwt_complex_double(x,y,n,cutoff) RESULT(is_comparable)
    COMPLEX(Double),           INTENT(IN) :: x, y
    INTEGER,                   INTENT(IN) :: n
    COMPLEX(Double), OPTIONAL, INTENT(IN) :: cutoff
    LOGICAL :: is_comparable
    COMPLEX(Double) :: c
    IF ( PRESENT(cutoff) ) THEN
      c = cutoff
    ELSE
      c = cmplx(dp_compare_cutoff,dp_compare_cutoff,double)
    END IF
    is_comparable = cwt_real_double(REAL(x,Double),REAL(y,Double),n,cutoff=real(c,double) ) .AND. &
                    cwt_real_double(aimag(x),aimag(y),n,cutoff=aimag(c))
  END FUNCTION cwt_complex_double

END MODULE compare_float_numbers