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

    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

    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))
    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

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 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