diff --git a/CHANGELOG.md b/CHANGELOG.md index da425d130..491a0e74e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,11 @@ Features available from the latest git source - new module `stdlib_io_npy` [#581](https://github.com/fortran-lang/stdlib/pull/581) - new procedures `save_npy`, `load_npy` +- update module `stdlib_math` + - new procedures `is_close` and `all_close` + [#488](https://github.com/fortran-lang/stdlib/pull/488) + - new procedures `arg`, `argd` and `argpi` + [#498](https://github.com/fortran-lang/stdlib/pull/498) Changes to existing modules diff --git a/doc/specs/stdlib_math.md b/doc/specs/stdlib_math.md index 2ce8bfce6..665d6ee5f 100644 --- a/doc/specs/stdlib_math.md +++ b/doc/specs/stdlib_math.md @@ -389,6 +389,129 @@ program demo_math_arange end program demo_math_arange ``` +### `arg` - Computes the phase angle in radian of a complex scalar + +#### Status + +Experimental + +#### Class + +Elemental function. + +#### Description + +`arg` computes the phase angle (radian version) of `complex` scalar in the interval (-π,π]. +The angles in `θ` are such that `z = abs(z)*exp((0.0, θ))`. + +#### Syntax + +`result = [[stdlib_math(module):arg(interface)]](z)` + +#### Arguments + +`z`: Shall be a `complex` scalar/array. +This is an `intent(in)` argument. + +#### Return value + +Returns the `real` type phase angle (radian version) of the `complex` argument `z`. + +Notes: Although the angle of the complex number `0` is undefined, `arg((0,0))` returns the value `0`. + +#### Example + +```fortran +program demo_math_arg + use stdlib_math, only: arg + print *, arg((0.0, 0.0)) !! 0.0 + print *, arg((3.0, 4.0)) !! 0.927 + print *, arg(2.0*exp((0.0, 0.5))) !! 0.5 +end program demo_math_arg +``` + +### `argd` - Computes the phase angle in degree of a complex scalar + +#### Status + +Experimental + +#### Class + +Elemental function. + +#### Description + +`argd` computes the phase angle (degree version) of `complex` scalar in the interval (-180.0,180.0]. +The angles in `θ` are such that `z = abs(z)*exp((0.0, θ*π/180.0))`. + +#### Syntax + +`result = [[stdlib_math(module):argd(interface)]](z)` + +#### Arguments + +`z`: Shall be a `complex` scalar/array. +This is an `intent(in)` argument. + +#### Return value + +Returns the `real` type phase angle (degree version) of the `complex` argument `z`. + +Notes: Although the angle of the complex number `0` is undefined, `argd((0,0))` returns the value `0`. + +#### Example + +```fortran +program demo_math_argd + use stdlib_math, only: argd + print *, argd((0.0, 0.0)) !! 0.0 + print *, argd((3.0, 4.0)) !! 53.1° + print *, argd(2.0*exp((0.0, 0.5))) !! 28.64° +end program demo_math_argd +``` + +### `argpi` - Computes the phase angle in circular of a complex scalar + +#### Status + +Experimental + +#### Class + +Elemental function. + +#### Description + +`argpi` computes the phase angle (IEEE circular version) of `complex` scalar in the interval (-1.0,1.0]. +The angles in `θ` are such that `z = abs(z)*exp((0.0, θ*π))`. + +#### Syntax + +`result = [[stdlib_math(module):argpi(interface)]](z)` + +#### Arguments + +`z`: Shall be a `complex` scalar/array. +This is an `intent(in)` argument. + +#### Return value + +Returns the `real` type phase angle (circular version) of the `complex` argument `z`. + +Notes: Although the angle of the complex number `0` is undefined, `argpi((0,0))` returns the value `0`. + +#### Example + +```fortran +program demo_math_argpi + use stdlib_math, only: argpi + print *, argpi((0.0, 0.0)) !! 0.0 + print *, argpi((3.0, 4.0)) !! 0.295 + print *, argpi(2.0*exp((0.0, 0.5))) !! 0.159 +end program demo_math_argpi +``` + ### `is_close` #### Description @@ -449,7 +572,6 @@ Returns a `logical` scalar/array. program demo_math_is_close use stdlib_math, only: is_close - use stdlib_error, only: check real :: x(2) = [1, 2], y, NAN y = -3 @@ -514,7 +636,6 @@ Returns a `logical` scalar. program demo_math_all_close use stdlib_math, only: all_close - use stdlib_error, only: check real :: x(2) = [1, 2], y, NAN complex :: z(4, 4) diff --git a/src/stdlib_math.fypp b/src/stdlib_math.fypp index 2db5d2a56..82f5961b2 100644 --- a/src/stdlib_math.fypp +++ b/src/stdlib_math.fypp @@ -14,7 +14,7 @@ module stdlib_math public :: EULERS_NUMBER_QP #:endif public :: DEFAULT_LINSPACE_LENGTH, DEFAULT_LOGSPACE_BASE, DEFAULT_LOGSPACE_LENGTH - public :: arange, is_close, all_close + public :: arange, arg, argd, argpi, is_close, all_close integer, parameter :: DEFAULT_LINSPACE_LENGTH = 100 integer, parameter :: DEFAULT_LOGSPACE_LENGTH = 50 @@ -27,6 +27,11 @@ module stdlib_math 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}$ @@ -296,6 +301,34 @@ module stdlib_math !> Version: experimental !> + !> `arg` computes the phase angle in the interval (-π,π]. + !> ([Specification](../page/specs/stdlib_math.html#arg)) + 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)) + 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)) + 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)) interface is_close @@ -341,6 +374,34 @@ contains #: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), z == (0.0_${k1}$, 0.0_${k1}$)) & + *180.0_${k1}$/PI_${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), z == (0.0_${k1}$, 0.0_${k1}$)) & + /PI_${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. @@ -361,4 +422,5 @@ contains end function gcd_${k1}$ #:endfor + end module stdlib_math diff --git a/src/tests/math/Makefile.manual b/src/tests/math/Makefile.manual index 9063a376b..5da73c5da 100644 --- a/src/tests/math/Makefile.manual +++ b/src/tests/math/Makefile.manual @@ -9,5 +9,4 @@ PROGS_SRC = test_linspace.f90 test_logspace.f90 \ $(SRCGEN): %.f90: %.fypp ../../common.fypp fypp -I../.. $(FYPPFLAGS) $< $@ - include ../Makefile.manual.test.mk diff --git a/src/tests/math/test_stdlib_math.fypp b/src/tests/math/test_stdlib_math.fypp index 72c3adf5b..4d8bb34a5 100644 --- a/src/tests/math/test_stdlib_math.fypp +++ b/src/tests/math/test_stdlib_math.fypp @@ -4,11 +4,15 @@ module test_stdlib_math use testdrive, only : new_unittest, unittest_type, error_type, check, skip_test - use stdlib_math, only: clip, is_close, all_close + use stdlib_math, only: clip, arg, argd, argpi, arange, is_close, all_close use stdlib_kinds, only: int8, int16, int32, int64, sp, dp, xdp, qp implicit none public :: collect_stdlib_math + + #:for k1 in REAL_KINDS + real(kind=${k1}$), parameter :: PI_${k1}$ = acos(-1.0_${k1}$) + #:endfor contains @@ -33,6 +37,13 @@ contains new_unittest("clip-real-quad", test_clip_rqp), & new_unittest("clip-real-quad-bounds", test_clip_rqp_bounds) & + !> Tests for arg/argd/argpi + #:for k1 in CMPLX_KINDS + , new_unittest("arg-cmplx-${k1}$", test_arg_${k1}$) & + , new_unittest("argd-cmplx-${k1}$", test_argd_${k1}$) & + , new_unittest("argpi-cmplx-${k1}$", test_argpi_${k1}$) & + #:endfor + !> Tests for `is_close` and `all_close` #:for k1 in REAL_KINDS , new_unittest("is_close-real-${k1}$", test_is_close_real_${k1}$) & @@ -211,7 +222,66 @@ contains #:endif end subroutine test_clip_rqp_bounds + + #:for k1 in CMPLX_KINDS + subroutine test_arg_${k1}$(error) + type(error_type), allocatable, intent(out) :: error + real(${k1}$), parameter :: tol = sqrt(epsilon(1.0_${k1}$)) + real(${k1}$), allocatable :: theta(:) + + #! For scalar + call check(error, abs(arg(2*exp((0.0_${k1}$, 0.5_${k1}$))) - 0.5_${k1}$) < tol, & + "test_nonzero_scalar") + if (allocated(error)) return + call check(error, abs(arg((0.0_${k1}$, 0.0_${k1}$)) - 0.0_${k1}$) < tol, & + "test_zero_scalar") + + #! and for array (180.0° see scalar version) + theta = arange(-179.0_${k1}$, 179.0_${k1}$, 3.58_${k1}$) + call check(error, all(abs(arg(exp(cmplx(0.0_${k1}$, theta/180*PI_${k1}$, ${k1}$))) - theta/180*PI_${k1}$) < tol), & + "test_array") + + end subroutine test_arg_${k1}$ + + subroutine test_argd_${k1}$(error) + type(error_type), allocatable, intent(out) :: error + real(${k1}$), parameter :: tol = sqrt(epsilon(1.0_${k1}$)) + real(${k1}$), allocatable :: theta(:) + + #! For scalar + call check(error, abs(argd((-1.0_${k1}$, 0.0_${k1}$)) - 180.0_${k1}$) < tol, & + "test_nonzero_scalar") + if (allocated(error)) return + call check(error, abs(argd((0.0_${k1}$, 0.0_${k1}$)) - 0.0_${k1}$) < tol, & + "test_zero_scalar") + + #! and for array (180.0° see scalar version) + theta = arange(-179.0_${k1}$, 179.0_${k1}$, 3.58_${k1}$) + call check(error, all(abs(argd(exp(cmplx(0.0_${k1}$, theta/180*PI_${k1}$, ${k1}$))) - theta) < tol), & + "test_array") + + end subroutine test_argd_${k1}$ + subroutine test_argpi_${k1}$(error) + type(error_type), allocatable, intent(out) :: error + real(${k1}$), parameter :: tol = sqrt(epsilon(1.0_${k1}$)) + real(${k1}$), allocatable :: theta(:) + + #! For scalar + call check(error, abs(argpi((-1.0_${k1}$, 0.0_${k1}$)) - 1.0_${k1}$) < tol, & + "test_nonzero_scalar") + if (allocated(error)) return + call check(error, abs(argpi((0.0_${k1}$, 0.0_${k1}$)) - 0.0_${k1}$) < tol, & + "test_zero_scalar") + + #! and for array (180.0° see scalar version) + theta = arange(-179.0_${k1}$, 179.0_${k1}$, 3.58_${k1}$) + call check(error, all(abs(argpi(exp(cmplx(0.0_${k1}$, theta/180*PI_${k1}$, ${k1}$))) - theta/180) < tol), & + "test_array") + + end subroutine test_argpi_${k1}$ + #:endfor + #:for k1 in REAL_KINDS subroutine test_is_close_real_${k1}$(error) type(error_type), allocatable, intent(out) :: error