[stdlib_math] add arg/argd/argpi

[zoziha]

• add arg/argd/argpi functions
• new style tests.

Description

(idea issue discussions see #489 )
arg computes the phase angle (radian version) of complex scalar in the interval (-π,π].
The angles in theta are such that z = abs(z)*exp((0.0, theta)).
arg API is:

result = [[stdlib_math(module):arg(interface)]](z)

Referring to the Fortran202X, I added argd/argpi.

(see 19-203r1.txt and 19-204r1.txt)

Prior Art

1. forlab/angle
2. matlab/angle
3. numpy/angle

[fiolj]

I don't know if there is an issue holding back this PR. I've looked into the files and I think they are fine.
The only thing that perhaps may be added is the inclusion of a use case with arrays in the tests. May be replacing in test_math.arg the line:

 print *, argpi(2.0*exp((0.0, PI))), argpi(2.0*exp(cmplx(0.0, -PI))) !! -0.999999940, 0.999999940

by something like:

 print *, argpi([2.0*exp((0.0, PI)), 2.0*exp(cmplx(0.0, -PI))]) !! -0.999999940, 0.999999940

[jvdp1]

I don't know if there is an issue holding back this PR.

This branch has some conflicts, but most likely a lack of reviewers too. I will try to review it in the next days.

[zoziha]

Thank you for your review @fiolj , I will use testdirve next time to refactor these test programs and solve these conficts.
Thank you for your review @jvdp1 , in fact, this PR content arg/angle is very easy to implement and review, before or after I modify the test programs, you can put forward your recommended changes.

[zoziha]

# Only 1 error: MacOS gfortran 9~11 real(sp) argd failed Starting argd-cmplx-sp ... (16/26) 179.999985 ... argd-cmplx-sp [FAILED] Message: test_nonzero_scalar

Very strange, gfortran on MacOS calculated the phase angle of single-precision cmplx(-1, 0)(only argd of real(sp) type on MacOS) and got 179.999985! It should be 180.0, but other platforms and other precisions are in line with expectations.

Should I increase the tolerance value, now tol=epsilon(..)?
Let tol=sqrt(epsilon(..))?

[fiolj]

# Only 1 error: MacOS gfortran 9~11 real(sp) argd failed Starting argd-cmplx-sp ... (16/26) 179.999985 ... argd-cmplx-sp [FAILED] Message: test_nonzero_scalar

Very strange, gfortran on MacOS calculated the phase angle of single-precision cmplx(-1, 0)(only argd of real(sp) type on MacOS) and got 179.999985! It should be 180.0, but other platforms and other precisions are in line with expectations.

Should I increase the tolerance value, now tol=epsilon(..)? Let tol=sqrt(epsilon(..))?

I can't see what is the problem, and cannot test on that operating system.
I am not sure that I understand correctly how it works, but if there is a difference (even if it were smaller than tol) and we scale it by 180/PI the difference will be larger than tol.
This, however, does not explain the observed difference in this case.

Being very simple I would try to see what values are giving in MacOs the expressions involved.
Could somebody try a program like the following? May be we can understand what is happening

!> test_argd program test_argd implicit none !> Single precision real numbers integer, parameter :: sp = selected_real_kind(6) real(sp), parameter :: tol = epsilon(1.0_sp) complex(sp) :: z = (-1.0_sp, 0.0_sp) real(sp) :: angle real(kind=sp), parameter :: PI = 4.0_sp*atan(1.0_sp) angle = aimag(log(z)) print *, 'PI =', PI print *, 'angle =', angle print *, 'Atan(z) =', atan(aimag(z), real(z)) print *, 'tol =', tol print *, 'angle=', angle*180/PI, 'diff=', 180._sp - angle*180/PI ! print *, 'obtained error =', abs(180._sp - 179.999985_sp) print *, 'original error =', abs(180._sp - 179.999985_sp)*PI/180._sp end program test_argd 

[jvdp1]

LGTM. I have only a major question regarding the tests: do you plan only to test for zero values?

EDIT:
My mistake: you also test for one non-zero value. Question: is it enough?

[fiolj]

I apologize for being so late with this comment, but there is a specific reason to use logarithms to calculate the angle? It seems unnecessary to evaluate the logarithm of the real part, if we are only keeping the angle.

If I am not missing something may be we can change the code to something like:

1 file changed, 3 insertions(+), 3 deletions(-) src/stdlib_math.fypp | 6 +++--- modified src/stdlib_math.fypp @@ -348,7 +348,7 @@ contains ${t1}$, intent(in) :: z real(${k1}$) :: result - result = aimag(log(z)) + result = atan2(aimag(z), real(z)) end function arg_${k1}$ @@ -356,7 +356,7 @@ contains ${t1}$, intent(in) :: z real(${k1}$) :: result - result = aimag(log(z))*180.0_${k1}$/PI_${k1}$ + result = atan2(aimag(z), real(z))*180.0_${k1}$/PI_${k1}$ end function argd_${k1}$ @@ -364,7 +364,7 @@ contains ${t1}$, intent(in) :: z real(${k1}$) :: result - result = aimag(log(z))/PI_${k1}$ + result = atan2(aimag(z), real(z))/PI_${k1}$ end function argpi_${k1}$ #:endfor 

[zoziha]

In fact, it is also possible to use atan2. I did not test the specific performance. If we use atan2, we need to add another judgment:

result = merge(0.0_rk, atan2(z%im, z%re), abs(z) == 0.0_rk) ! or result = merge(0.0_rk, atan2(z%im, z%re), z%im == 0.0_rk .and. z%re == 0.0_rk)

print *, log(cmplx(0.0,0.0)) !! (-Infinity,0.0000000E+00) (IFORT compiler)

[zoziha]

There may indeed be an unstable place here:
It can be found in the compilation check of GFortran:

app/example.f90:8:17: 8 | print *, log(cmplx(0.0,0.0)) | 1 Error: Complex argument of LOG at (1) cannot be zero

But at runtime, log(cmplx(0.0,0.0)) can be calculated. Perhaps this is one of the reasons for using atan2.

The main reason I use aimag(log(..)) is that it is more concise and intuitive, while numpy/angle uses arctan2(..) (like Fortran's atan2) (see numpy/angle API, and numpy/angle source).

[zoziha]

I tried to use the atan2 solution, but CI of macos still failed. It seems that it might be a bug in gfortran on macos?

[gareth-nx]

Thanks for this -- a few changes suggested here -- unrelated to the macos issue.

[zoziha]

Three main changes:

• Internal implementation: aimag(log(..)) -> atan2(..)
• Tolerance of unit tests: epsilon(..) -> sqrt(epsilon(..)) (Since gfortran on MacOS failed, maybe a bug, see [stdlib_math] add arg/argd/argpi #498 (comment))
• Updated testing of the array content for more angles.
• ...

[gareth-nx]

Looks good to me, thanks.

[fiolj]

* Tolerance of unit tests: `epsilon(..) -> sqrt(epsilon(..))` (Since gfortran on MacOS failed, maybe a bug) 

Changes look good. Minor issue: using sqrt(epsilon) for tolerance is probably too lax, especially for double precision where typical errors (for example in the values of the tests) are aprox 1.e-14 and tolerance is about 1.e-8. Probably something like
tol = A*epsilon(...) with values of A=500, or A=1000, would be more adequate.
Anyway, as it is looks good and, in my opinion, could be merged

[awvwgk]

Looks good to me. Do you mind adding this addition to the CHANGELOG.md file?

[zoziha]

Fine, I will add later.

[zoziha]

CHANGELOG.md has been updated, and an unnecessary addition of **_close examples in stdlib_math.md has been removed.

[awvwgk]

Thanks a lot, will merge once the CI has passed.