# basic_float_math.F90 Source File

## Source Code

#include "macros.h"

module basic_float_math
use constants, only : sp, dp, EPS
implicit none

private
public :: isclose, operator(.isclose.), near_zero, conjgt, get_nan, is_nan

interface operator(.isclose.)
module procedure isclose_for_operator_real_dp
module procedure isclose_for_operator_real_sp
module procedure isclose_for_operator_complex_dp
module procedure isclose_for_operator_complex_sp
end interface

interface isclose
module procedure isclose_real_dp
module procedure isclose_real_sp
module procedure isclose_complex_dp
module procedure isclose_complex_sp
end interface

interface near_zero
module procedure near_zero_real_dp
module procedure near_zero_real_sp
module procedure near_zero_complex_dp
module procedure near_zero_complex_sp
end interface

interface conjgt
module procedure conjgt_real_dp
module procedure conjgt_real_sp
module procedure conjgt_complex_dp
module procedure conjgt_complex_sp
end interface

contains

!> @brief
!> Compare floating point numbers for equality
!>
!> @details
!> The comparison is symmetric and decorrelates atol and rtol.
!> A good discussion can be found here
!> https://github.com/numpy/numpy/issues/10161
!>
!> @param[in] a
!> @param[in] b
!> @param[in] atol The absolute tolerance. Defaults to 0.
!> @param[in] rtol The relative tolerance. Defaults to 1e-9.
elemental function isclose_real_dp(a, b, atol, rtol) result(res)
real(dp), intent(in) :: a, b
real(dp), intent(in), optional :: atol, rtol
logical :: res

real(dp) :: atol_, rtol_

def_default(rtol_, rtol, 1e-9_dp)
def_default(atol_, atol, 0._dp)

res = abs(a - b) <= max(atol_, rtol_ * min(abs(a), abs(b)))
end function

!> Operator functions may only have two arguments.
elemental function isclose_for_operator_real_dp(a, b) result(res)
real(dp), intent(in) :: a, b
logical :: res

res = isclose(a, b, atol=EPS)
end function

elemental function near_zero_real_dp(x, epsilon) result(res)
real(dp), intent(in) :: x
real(dp), intent(in), optional :: epsilon
logical :: res

real(dp) :: epsilon_

def_default(epsilon_, epsilon, EPS)

res = isclose(x, 0._dp, atol=epsilon_, rtol=0._dp)
end function

!> @brief
!> Compare floating point numbers for equality
!>
!> @details
!> The comparison is symmetric and decorrelates atol and rtol.
!> A good discussion can be found here
!> https://github.com/numpy/numpy/issues/10161
!>
!> @param[in] a
!> @param[in] b
!> @param[in] atol The absolute tolerance. Defaults to 0.
!> @param[in] rtol The relative tolerance. Defaults to 1e-9.
elemental function isclose_real_sp(a, b, atol, rtol) result(res)
real(sp), intent(in) :: a, b
real(dp), intent(in), optional :: atol, rtol
logical :: res

real(dp) :: atol_, rtol_

def_default(rtol_, rtol, 1e-9_dp)
def_default(atol_, atol, 0._dp)

res = abs(a - b) <= max(atol_, rtol_ * min(abs(a), abs(b)))
end function

!> Operator functions may only have two arguments.
elemental function isclose_for_operator_real_sp(a, b) result(res)
real(sp), intent(in) :: a, b
logical :: res

res = isclose(a, b, atol=EPS)
end function

elemental function near_zero_real_sp(x, epsilon) result(res)
real(sp), intent(in) :: x
real(dp), intent(in), optional :: epsilon
logical :: res

real(dp) :: epsilon_

def_default(epsilon_, epsilon, EPS)

res = isclose(x, 0._sp, atol=epsilon_, rtol=0._dp)
end function

!> @brief
!> Compare floating point numbers for equality
!>
!> @details
!> The comparison is symmetric and decorrelates atol and rtol.
!> A good discussion can be found here
!> https://github.com/numpy/numpy/issues/10161
!>
!> @param[in] a
!> @param[in] b
!> @param[in] atol The absolute tolerance. Defaults to 0.
!> @param[in] rtol The relative tolerance. Defaults to 1e-9.
elemental function isclose_complex_dp(a, b, atol, rtol) result(res)
complex(dp), intent(in) :: a, b
real(dp), intent(in), optional :: atol, rtol
logical :: res

real(dp) :: atol_, rtol_

def_default(rtol_, rtol, 1e-9_dp)
def_default(atol_, atol, 0._dp)

res = abs(a - b) <= max(atol_, rtol_ * min(abs(a), abs(b)))
end function

!> Operator functions may only have two arguments.
elemental function isclose_for_operator_complex_dp(a, b) result(res)
complex(dp), intent(in) :: a, b
logical :: res

res = isclose(a, b, atol=EPS)
end function

elemental function near_zero_complex_dp(x, epsilon) result(res)
complex(dp), intent(in) :: x
real(dp), intent(in), optional :: epsilon
logical :: res

real(dp) :: epsilon_

def_default(epsilon_, epsilon, EPS)

res = isclose(x, (0._dp, 0._dp), atol=epsilon_, rtol=0._dp)
end function

!> @brief
!> Compare floating point numbers for equality
!>
!> @details
!> The comparison is symmetric and decorrelates atol and rtol.
!> A good discussion can be found here
!> https://github.com/numpy/numpy/issues/10161
!>
!> @param[in] a
!> @param[in] b
!> @param[in] atol The absolute tolerance. Defaults to 0.
!> @param[in] rtol The relative tolerance. Defaults to 1e-9.
elemental function isclose_complex_sp(a, b, atol, rtol) result(res)
complex(sp), intent(in) :: a, b
real(dp), intent(in), optional :: atol, rtol
logical :: res

real(dp) :: atol_, rtol_

def_default(rtol_, rtol, 1e-9_dp)
def_default(atol_, atol, 0._dp)

res = abs(a - b) <= max(atol_, rtol_ * min(abs(a), abs(b)))
end function

!> Operator functions may only have two arguments.
elemental function isclose_for_operator_complex_sp(a, b) result(res)
complex(sp), intent(in) :: a, b
logical :: res

res = isclose(a, b, atol=EPS)
end function

elemental function near_zero_complex_sp(x, epsilon) result(res)
complex(sp), intent(in) :: x
real(dp), intent(in), optional :: epsilon
logical :: res

real(dp) :: epsilon_

def_default(epsilon_, epsilon, EPS)

res = isclose(x, (0._sp, 0._sp), atol=epsilon_, rtol=0._dp)
end function

elemental function conjgt_real_dp(x) result(res)
real(dp), intent(in) :: x
real(dp) :: res
res = x
end function
elemental function conjgt_real_sp(x) result(res)
real(sp), intent(in) :: x
real(sp) :: res
res = x
end function

elemental function conjgt_complex_dp(x) result(res)
complex(dp), intent(in) :: x
complex(dp) :: res
res = conjg(x)
end function
elemental function conjgt_complex_sp(x) result(res)
complex(sp), intent(in) :: x
complex(sp) :: res
res = conjg(x)
end function

real(dp) pure function get_nan()
! If all of the compilers supported ieee_arithmetic
! --> could use ieee_value(1.0_dp, ieee_quiet_nan)
real(dp) :: a, b
a = 1
b = 1
get_nan = log(a - 2 * b)
end function

elemental logical function is_nan(r)
! If all of the compilers supported ieee_arithmetic
! --> could use ieee_is_nan (r)
real(dp), intent(in) :: r

#ifdef GFORTRAN_
is_nan = isnan(r)
#else
is_nan = r /= r
#endif
end function

end module