#: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 :: stdlib_meshgrid_ij, stdlib_meshgrid_xy
public :: arange, arg, argd, argpi, is_close, all_close, diff, meshgrid
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
!> Values for optional argument `indexing` of `meshgrid`
integer, parameter :: stdlib_meshgrid_xy = 0, stdlib_meshgrid_ij = 1
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
!> Version: experimental
!>
!> Computes a list of coordinate matrices from coordinate vectors.
!> ([Specification](../page/specs/stdlib_math.html#meshgrid))
interface meshgrid
#:set RANKS = range(1, MAXRANK + 1)
#:for k1, t1 in IR_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("meshgrid", rank, t1, k1)
module subroutine ${RName}$(&
${"".join(f"x{i}, " for i in range(1, rank + 1))}$ &
${"".join(f"xm{i}, " for i in range(1, rank + 1))}$ &
indexing &
)
#:for i in range(1, rank + 1)
${t1}$, intent(in) :: x${i}$(:)
${t1}$, intent(out) :: xm${i}$ ${ranksuffix(rank)}$
#:endfor
integer, intent(in), optional :: indexing
end subroutine ${RName}$
#:endfor
#:endfor
end interface meshgrid
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
