stdlib_math.fypp

#:include "common.fypp"
#:set IR_KINDS_TYPES = INT_KINDS_TYPES + REAL_KINDS_TYPES
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES

module stdlib_math
    use stdlib_kinds, only: int8, int16, int32, int64, sp, dp, xdp, qp
    use stdlib_optval, only: optval

    implicit none
    private
    public :: clip, gcd, linspace, logspace
    public :: EULERS_NUMBER_SP, EULERS_NUMBER_DP
#:if WITH_QP
    public :: EULERS_NUMBER_QP
#:endif
    public :: DEFAULT_LINSPACE_LENGTH, DEFAULT_LOGSPACE_BASE, DEFAULT_LOGSPACE_LENGTH
    public :: arange, arg, argd, argpi, is_close, all_close, diff

    integer, parameter :: DEFAULT_LINSPACE_LENGTH = 100
    integer, parameter :: DEFAULT_LOGSPACE_LENGTH = 50
    integer, parameter :: DEFAULT_LOGSPACE_BASE = 10

    ! Useful constants for lnspace
    real(sp), parameter :: EULERS_NUMBER_SP = exp(1.0_sp)
    real(dp), parameter :: EULERS_NUMBER_DP = exp(1.0_dp)
#:if WITH_QP
    real(qp), parameter :: EULERS_NUMBER_QP = exp(1.0_qp)
#:endif

    !> Useful constants `PI` for `argd/argpi`
    #:for k1 in REAL_KINDS
    real(kind=${k1}$), parameter :: PI_${k1}$ = acos(-1.0_${k1}$)
    #:endfor

    interface clip
    #:for k1, t1 in IR_KINDS_TYPES
        module procedure clip_${k1}$
    #:endfor
    end interface clip

    !> Returns the greatest common divisor of two integers
    !> ([Specification](../page/specs/stdlib_math.html#gcd))
    !>
    !> Version: experimental
    interface gcd
    #:for k1, t1 in INT_KINDS_TYPES
        module procedure gcd_${k1}$
    #:endfor
    end interface gcd

    interface linspace
    !! Version: Experimental
    !!
    !! Create rank 1 array of linearly spaced elements
    !! If the number of elements is not specified, create an array with size 100. If n is a negative value,
    !! return an array with size 0. If n = 1, return an array whose only element is end
    !!([Specification](../page/specs/stdlib_math.html#linspace-create-a-linearly-spaced-rank-one-array))
    #:for k1, t1 in RC_KINDS_TYPES
      #:set RName = rname("linspace_default", 1, t1, k1)
      pure module function ${RName}$(start, end) result(res)
        ${t1}$, intent(in) :: start
        ${t1}$, intent(in) :: end

        ${t1}$ :: res(DEFAULT_LINSPACE_LENGTH)
      end function ${RName}$
    #:endfor

    #:for k1, t1 in RC_KINDS_TYPES
      #:set RName = rname("linspace_n", 1, t1, k1)
      pure module function ${RName}$(start, end, n) result(res)
        ${t1}$, intent(in) :: start
        ${t1}$, intent(in) :: end
        integer, intent(in) :: n

        ${t1}$ :: res(max(n, 0))
      end function ${RName}$
    #:endfor


    ! Add support for integer linspace
    !!
    !! When dealing with integers as the `start` and `end` parameters, the return type is always a `real(dp)`.
    #:for k1, t1 in INT_KINDS_TYPES
      #:set RName = rname("linspace_default", 1, t1, k1)
      #! The interface for INT_KINDS_TYPES cannot be combined with RC_KINDS_TYPES
      #! because the output for integer types is always a real with dp.
      pure module function ${RName}$(start, end) result(res)
        ${t1}$, intent(in) :: start
        ${t1}$, intent(in) :: end

        real(dp) :: res(DEFAULT_LINSPACE_LENGTH)
      end function ${RName}$
    #:endfor

    #:for k1, t1 in INT_KINDS_TYPES
      #:set RName = rname("linspace_n", 1, t1, k1)
      pure module function ${RName}$(start, end, n) result(res)
        ${t1}$, intent(in) :: start
        ${t1}$, intent(in) :: end
        integer, intent(in) :: n

        real(dp) :: res(max(n, 0))
      end function ${RName}$
    #:endfor

  end interface

  interface logspace
  !! Version: Experimental
  !!
  !! Create rank 1 array of logarithmically spaced elements from base**start to base**end.
  !! If the number of elements is not specified, create an array with size 50. If n is a negative value,
  !! return an array with size 0. If n = 1, return an array whose only element is base**end. If no base
  !! is specified, logspace will default to using a base of 10
  !!
  !!([Specification](../page/specs/stdlib_math.html#logspace-create-a-logarithmically-spaced-rank-one-array))
  #!=========================================================
  #!=     logspace(start, end)                              =
  #!=========================================================
  #:for k1, t1 in RC_KINDS_TYPES
    #:set RName = rname("logspace", 1, t1, k1, "default")
    pure module function ${RName}$(start, end) result(res)

      ${t1}$, intent(in) :: start
      ${t1}$, intent(in) :: end

      ${t1}$ :: res(DEFAULT_LOGSPACE_LENGTH)

    end function ${RName}$
  #:endfor
  #! Integer support
  #:set RName = rname("logspace", 1, "integer(int32)", "int32", "default")
    pure module function ${RName}$(start, end) result(res)

      integer, intent(in) :: start
      integer, intent(in) :: end

      real(dp) :: res(DEFAULT_LOGSPACE_LENGTH)

  end function ${RName}$

  #!=========================================================
  #!=     logspace(start, end, n)                           =
  #!=========================================================
  #:for k1, t1 in RC_KINDS_TYPES
    #:set RName = rname("logspace", 1, t1, k1, "n")
    pure module function ${RName}$(start, end, n) result(res)
      ${t1}$, intent(in) :: start
      ${t1}$, intent(in) :: end
      integer, intent(in) :: n

      ${t1}$ :: res(max(n, 0))
    end function ${RName}$
  #:endfor
  #! Integer support
  #:set RName = rname("logspace", 1, "integer(int32)", "int32", "n")
    pure module function ${RName}$(start, end, n) result(res)
      integer, intent(in) :: start
      integer, intent(in) :: end
      integer, intent(in) :: n

      real(dp) :: res(n)
    end function ${RName}$

  #!=========================================================
  #!=     logspace(start, end, n, base)                     =
  #!=========================================================
  #! Need another function where base is not optional,
  #! otherwise the compiler can not differentiate between
  #! generic calls to logspace_n where a base is not present
  #! ========================================================
  #:for k1, t1 in REAL_KINDS_TYPES
    ! Generate logarithmically spaced sequence from ${k1}$ base to the powers
    ! of ${k1}$ start and end. [base^start, ... , base^end]
    ! Different combinations of parameter types will lead to different result types.
    ! Those combinations are indicated in the body of each function.
    #:set RName = rname("logspace", 1, t1, k1, "n_rbase")
    pure module function ${RName}$(start, end, n, base) result(res)
      ${t1}$, intent(in) :: start
      ${t1}$, intent(in) :: end
      integer, intent(in) :: n
      ${t1}$, intent(in) :: base
      ! real(${k1}$) endpoints + real(${k1}$) base = real(${k1}$) result
      ${t1}$ :: res(max(n, 0))
    end function ${RName}$

    #:set RName = rname("logspace", 1, t1, k1, "n_cbase")
    pure module function ${RName}$(start, end, n, base) result(res)
      ${t1}$, intent(in) :: start
      ${t1}$, intent(in) :: end
      integer, intent(in) :: n
      complex(${k1}$), intent(in) :: base
      ! real(${k1}$) endpoints + complex(${k1}$) base = complex(${k1}$) result
      ${t1}$ :: res(max(n, 0))
    end function ${RName}$

    #:set RName = rname("logspace", 1, t1, k1, "n_ibase")
    pure module function ${RName}$(start, end, n, base) result(res)
      ${t1}$, intent(in) :: start
      ${t1}$, intent(in) :: end
      integer, intent(in) :: n
      integer, intent(in) :: base
      ! real(${k1}$) endpoints + integer base = real(${k1}$) result
      ${t1}$ :: res(max(n, 0))
    end function ${RName}$
  #:endfor
  #! ========================================================
  #! ========================================================
  #:for k1, t1 in CMPLX_KINDS_TYPES
    ! Generate logarithmically spaced sequence from ${k1}$ base to the powers
    ! of ${k1}$ start and end. [base^start, ... , base^end]
    ! Different combinations of parameter types will lead to different result types.
    ! Those combinations are indicated in the body of each function.
    #:set RName = rname("logspace", 1, t1, k1, "n_rbase")
    pure module function ${RName}$(start, end, n, base) result(res)
      ${t1}$, intent(in) :: start
      ${t1}$, intent(in) :: end
      integer, intent(in) :: n
      real(${k1}$), intent(in) :: base
      ! complex(${k1}$) endpoints + real(${k1}$) base = complex(${k1}$) result
      ${t1}$ :: res(max(n, 0))
    end function ${RName}$

    #:set RName = rname("logspace", 1, t1, k1, "n_cbase")
    pure module function ${RName}$(start, end, n, base) result(res)
      ${t1}$, intent(in) :: start
      ${t1}$, intent(in) :: end
      integer, intent(in) :: n
      complex(${k1}$), intent(in) :: base
      ! complex(${k1}$) endpoints + complex(${k1}$) base = complex(${k1}$) result
      ${t1}$ :: res(max(n, 0))
    end function ${RName}$

    #:set RName = rname("logspace", 1, t1, k1, "n_ibase")
    pure module function ${RName}$(start, end, n, base) result(res)
      ${t1}$, intent(in) :: start
      ${t1}$, intent(in) :: end
      integer, intent(in) :: n
      integer, intent(in) :: base
      ! complex(${k1}$) endpoints + integer base = complex(${k1}$) result
      ${t1}$ :: res(max(n, 0))
    end function ${RName}$
  #:endfor
  #! ========================================================
  #! ========================================================
  #! Provide support for Integer start/endpoints
    ! Generate logarithmically spaced sequence from ${k1}$ base to the powers
    ! of ${k1}$ start and end. [base^start, ... , base^end]
    ! Different combinations of parameter types will lead to different result types.
    ! Those combinations are indicated in the body of each function.
  #:for k1 in REAL_KINDS
    #:set RName = rname("logspace", 1, "integer(int32)", "int32", "n_r" + str(k1) + "base")
    pure module function ${RName}$(start, end, n, base) result(res)
      integer, intent(in) :: start
      integer, intent(in) :: end
      integer, intent(in) :: n
      real(${k1}$), intent(in) :: base
      ! integer endpoints + real(${k1}$) base = real(${k1}$) result
      real(${k1}$) :: res(max(n, 0))
    end function ${RName}$

    #:set RName = rname("logspace", 1, "integer(int32)", "int32", "n_c" + str(k1) + "base")
    pure module function ${RName}$(start, end, n, base) result(res)
      integer, intent(in) :: start
      integer, intent(in) :: end
      integer, intent(in) :: n
      complex(${k1}$), intent(in) :: base
      ! integer endpoints + complex(${k1}$) base = complex(${k1}$) result
      complex(${k1}$) :: res(max(n, 0))
    end function ${RName}$
  #:endfor

    #:set RName = rname("logspace", 1, "integer(int32)", "int32", "n_ibase")
    pure module function ${RName}$(start, end, n, base) result(res)
      integer, intent(in) :: start
      integer, intent(in) :: end
      integer, intent(in) :: n
      integer, intent(in) :: base
      ! integer endpoints + integer base = integer result
      integer :: res(max(n, 0))
    end function ${RName}$


  end interface

    !> Version: experimental
    !>
    !> `arange` creates a one-dimensional `array` of the `integer/real` type 
    !>  with fixed-spaced values of given spacing, within a given interval.
    !> ([Specification](../page/specs/stdlib_math.html#arange-function))
    interface arange
        #:set RI_KINDS_TYPES = REAL_KINDS_TYPES + INT_KINDS_TYPES
        #:for k1, t1 in RI_KINDS_TYPES
        pure module function arange_${t1[0]}$_${k1}$(start, end, step) result(result)
            ${t1}$, intent(in) :: start
            ${t1}$, intent(in), optional :: end, step
            ${t1}$, allocatable :: result(:)
        end function arange_${t1[0]}$_${k1}$
        #:endfor
    end interface arange

    !> Version: experimental
    !>
    !> `arg` computes the phase angle in the interval (-π,π].
    !> ([Specification](../page/specs/stdlib_math.html#arg-function))
    interface arg
        #:for k1 in CMPLX_KINDS
        procedure :: arg_${k1}$
        #:endfor
    end interface arg

    !> Version: experimental
    !>
    !> `argd` computes the phase angle of degree version in the interval (-180.0,180.0].
    !> ([Specification](../page/specs/stdlib_math.html#argd-function))
    interface argd
        #:for k1 in CMPLX_KINDS
        procedure :: argd_${k1}$
        #:endfor
    end interface argd

    !> Version: experimental
    !>
    !> `argpi` computes the phase angle of circular version in the interval (-1.0,1.0].
    !> ([Specification](../page/specs/stdlib_math.html#argpi-function))
    interface argpi
        #:for k1 in CMPLX_KINDS
        procedure :: argpi_${k1}$
        #:endfor
    end interface argpi
    
    !> Returns a boolean scalar/array where two scalar/arrays are element-wise equal within a tolerance.
    !> ([Specification](../page/specs/stdlib_math.html#is_close-function))
    interface is_close
        #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
        #:for k1, t1 in RC_KINDS_TYPES
        elemental module logical function is_close_${t1[0]}$${k1}$(a, b, rel_tol, abs_tol, equal_nan) result(close)
            ${t1}$, intent(in) :: a, b
            real(${k1}$), intent(in), optional :: rel_tol, abs_tol
            logical, intent(in), optional :: equal_nan
        end function is_close_${t1[0]}$${k1}$
        #:endfor
    end interface is_close

    !> Version: experimental
    !>
    !> Returns a boolean scalar where two arrays are element-wise equal within a tolerance.
    !> ([Specification](../page/specs/stdlib_math.html#all_close-function))
    interface all_close
        #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
        #:set RANKS = range(1, MAXRANK + 1)
        #:for k1, t1 in RC_KINDS_TYPES
        #:for r1 in RANKS
        logical pure module function all_close_${r1}$_${t1[0]}$${k1}$(a, b, rel_tol, abs_tol, equal_nan) result(close)
            ${t1}$, intent(in) :: a${ranksuffix(r1)}$, b${ranksuffix(r1)}$
            real(${k1}$), intent(in), optional :: rel_tol, abs_tol
            logical, intent(in), optional :: equal_nan
        end function all_close_${r1}$_${t1[0]}$${k1}$
        #:endfor
        #:endfor
    end interface all_close
    
    !> Version: experimental
    !>
    !> Computes differences between adjacent elements of an array.
    !> ([Specification](../page/specs/stdlib_math.html#diff-function))
    interface diff
        #:set RI_KINDS_TYPES = REAL_KINDS_TYPES + INT_KINDS_TYPES
        #:for k1, t1 in RI_KINDS_TYPES
        pure module function diff_1_${k1}$(x, n, prepend, append) result(y)
            ${t1}$, intent(in) :: x(:)
            integer, intent(in), optional :: n
            ${t1}$, intent(in), optional :: prepend(:), append(:)
            ${t1}$, allocatable :: y(:)
        end function diff_1_${k1}$
        pure module function diff_2_${k1}$(X, n, dim, prepend, append) result(y)
            ${t1}$, intent(in) :: x(:, :)
            integer, intent(in), optional :: n, dim
            ${t1}$, intent(in), optional :: prepend(:, :), append(:, :)
            ${t1}$, allocatable :: y(:, :)
        end function diff_2_${k1}$
        #:endfor
    end interface diff

contains

    #:for k1, t1 in IR_KINDS_TYPES
    elemental function clip_${k1}$(x, xmin, xmax) result(res)
        ${t1}$, intent(in) :: x
        ${t1}$, intent(in) :: xmin
        ${t1}$, intent(in) :: xmax
        ${t1}$ :: res

        res = max(min(x, xmax), xmin)
    end function clip_${k1}$

    #:endfor

    #:for k1, t1 in CMPLX_KINDS_TYPES
    elemental function arg_${k1}$(z) result(result) 
        ${t1}$, intent(in) :: z
        real(${k1}$) :: result

        result = merge(0.0_${k1}$, atan2(z%im, z%re), z == (0.0_${k1}$, 0.0_${k1}$))

    end function arg_${k1}$

    elemental function argd_${k1}$(z) result(result) 
        ${t1}$, intent(in) :: z
        real(${k1}$) :: result

        result = merge(0.0_${k1}$, atan2(z%im, z%re)*180.0_${k1}$/PI_${k1}$, &
                 z == (0.0_${k1}$, 0.0_${k1}$))

    end function argd_${k1}$

    elemental function argpi_${k1}$(z) result(result) 
        ${t1}$, intent(in) :: z
        real(${k1}$) :: result

        result = merge(0.0_${k1}$, atan2(z%im, z%re)/PI_${k1}$, &
                 z == (0.0_${k1}$, 0.0_${k1}$))
                 

    end function argpi_${k1}$
    #:endfor

    #:for k1, t1 in INT_KINDS_TYPES
    !> Returns the greatest common divisor of two integers of kind ${k1}$
    !> using the Euclidean algorithm.
    elemental function gcd_${k1}$(a, b) result(res)
        ${t1}$, intent(in) :: a
        ${t1}$, intent(in) :: b
        ${t1}$ :: res

        ${t1}$ :: rem, tmp

        rem = min(abs(a), abs(b))
        res = max(abs(a), abs(b))
        do while (rem /= 0_${k1}$)
          tmp = rem
          rem = mod(res, rem)
          res = tmp
        end do
    end function gcd_${k1}$

    #:endfor
    
end module stdlib_math