[stdlib_math] add arg/argd/argpi (#498)

Author: zoziha

CHANGELOG.md

@@ -16,6 +16,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
 

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)
     

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

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

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