basic_float_math.F90 Source File

#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