 |
forums.silverfrost.com
Welcome to the Silverfrost forums
|
| View previous topic :: View next topic |
| Author |
Message |
Kenneth_Smith
Joined: 18 May 2012
Posts: 797
Location: Hamilton, Lanarkshire, Scotland.
|
 Posted: Thu Feb 27, 2025 10:53 am Post subject: Puzzle asin(z) and acos(z) |
 |
|
The complex arcsin and arccos functions can each be expressed in two ways:
| Code: | arcsin(z) = -j*log( sqrt(1.0 - z*z) + j*z ) Eqn 1
= j*log( sqrt(1.0 - z*z) - j*z ) Eqn 2
arccos(z) = -j*log( j*sqrt(1.0 � z*z) + z ) Eqn 3
= pi/2 � arcsin(z) Eqn 4
|
The square root term 1.0 � z*z is common to three of these, and the addition of +j*z or -j*z in eqns 1 and 2 determine what is passed to the LOG function for the arcsin.
Under certain conditions this could result in both the real and imaginary parts of the argument to the log function being zero.
Now consider the following program:
| Code: | program p
use iso_fortran_env
implicit none
integer, parameter :: wp = kind(1.d0)
complex(kind=wp) :: z1,z2
print*, compiler_version()
print*
z1 = cmplx(500.0D6,500.0D6,kind=wp)
z2 = cmplx(500.0D6,-500.0D6,kind=wp)
print*, acos(z1)
print*, acos(z2) !Real and imaginary parts of argument to LOG are both zero
print*, asin(z1) !Real and imaginary parts of argument to LOG are both zero
print*, asin(z2)
end program p
|
With the problematic lines commented out, trying three different compilers:
| Code: | FTN95 v9.10.0
(0.78539816339,-21.0698394272)
(0.78539816339,-21.0698394272)
GCC version 12.3.0
(1.5707963267948966,-Inf)
(0.78539816339744828,-21.069839427226384)
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel
(R) 64, Version 2021.4.0 Build 20210910_000000
(0.785398163397448,-21.0698394272264)
(0.785398163397448,-21.0698394272264)
|
We see that FTN95 returns the same values as Ifort, but gFortran does something different for ACOS(Z1).
Last edited by Kenneth_Smith on Thu Feb 27, 2025 11:00 am; edited 1 time in total |
|
| Back to top |
|
 |
Kenneth_Smith
Joined: 18 May 2012
Posts: 797
Location: Hamilton, Lanarkshire, Scotland.
|
| Posted: Thu Feb 27, 2025 10:57 am Post subject: |
|
|
Continued .....
For the full code above, with the lines problematic to FTN not commented out, we get:
| Code: | GCC version 12.3.0
(1.5707963267948966,-Inf)
(0.78539816339744828,21.069839427226384)
(0.0000000000000000,Inf)
(0.78539816339744828,-21.069839427226384)
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel
(R) 64, Version 2021.4.0 Build 20210910_000000
(0.785398163397448,-21.0698394272264)
(0.785398163397448,21.0698394272264)
(0.785398163397448,21.0698394272264)
(0.785398163397448,-21.0698394272264)
|
I think FTN95 uses just one of the two available functions for asin(z) and acos(z), while Ifort selects which equation to use based on the sign of the imaginary part of Z.
Using the sign of the imaginary part of Z may not be sufficient, since there are issues with logarithm branch cuts (when the complex number passed to the log function is a negative real number). If the imaginary part of z is very large the imaginary part goes towards + infinity, and � minus infinity for very large negative imaginary part. This extra complication is being considered by gFortran.
It seems that the performance of FTN95 could be improved here, or a user defined function defined for these functions, but what should be the benchmark to compare with?
I�ve pushed the envelope a bit here in terms of my knowledge of maths, but perhaps other forum reader can see something that I have missed. |
|
| Back to top |
|
|
PaulLaidler
Site Admin
Joined: 21 Feb 2005
Posts: 8165
Location: Salford, UK
|
| Posted: Thu Feb 27, 2025 12:41 pm Post subject: |
|
|
Ken
I have made a note to look into this. |
|
| Back to top |
|
|
Kenneth_Smith
Joined: 18 May 2012
Posts: 797
Location: Hamilton, Lanarkshire, Scotland.
|
| Posted: Thu Feb 27, 2025 1:58 pm Post subject: |
|
|
No worries Paul, this is not urgent - more of a frustrating curiosity.
The single precision SLATEC routines CASIN(z) and CACOS(z) align with iFort for z = cmplx(500.0e6,500.0e6), but differ for z = cmplx(500.0e6,-500.0e6)
| Code: | Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel
(R) 64, Version 2021.4.0 Build 20210910_000000
Test asin of (5.0000000E+08,5.0000000E+08)
From SLATEC (0.7853982,21.06984)
From compiler (0.7853982,21.06984)
Test acos of (5.0000000E+08,5.0000000E+08)
From SLATEC (0.7853981,-21.06984)
From compiler (0.7853982,-21.06984)
Test asin of (5.0000000E+08,-5.0000000E+08)
From SLATEC (-1.570796,3.465736)
From compiler (0.7853982,-21.06984)
Test acos of (5.0000000E+08,-5.0000000E+08)
From SLATEC (3.141593,-3.465736)
From compiler (0.7853982,21.06984)
|
There is no alignment between SLATEC and gFortran
| Code: | GCC version 12.3.0
Test asin of (500000000.,500000000.)
From SLATEC (0.785398245,21.0698395)
From compiler (0.00000000,Inf)
Test acos of (500000000.,500000000.)
From SLATEC (0.785398126,-21.0698395)
From compiler (1.57079637,-Inf)
Test asin of (500000000.,-500000000.)
From SLATEC (-1.57079625,3.46573591)
From compiler (0.785398185,-21.0698395)
Test acos of (500000000.,-500000000.)
From SLATEC (3.14159250,-3.46573591)
From compiler (0.785398185,21.0698395
|
Rather a big can of worms! |
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 146
Location: United Kingdom
|
| Posted: Sat Mar 01, 2025 3:38 pm Post subject: |
|
|
Hi Ken
I tried a test code and get some values close to the expected result but it falls over for large values as a complex input to acos.
Using the explicit functions for acos and asin as functions:
| Code: |
module complex_acos_asin
implicit none
contains
function cacos_f(z) result(c_acos)
implicit none
complex*8, intent(in) :: z
complex*8 :: c_acos, unit_i=(0.0d0, 1.0d0)
! cacos(z) = -i * log(z + i * sqrt(1 - z^2))
c_acos = - unit_i * log(z + unit_i * sqrt(1.0d0 - z**2))
end function cacos_f
function casin_f(z) result(c_asin)
implicit none
complex*8, intent(in) :: z
complex*8 :: c_asin, unit_i=(0.0d0, 1.0d0)
! casin(z) = -i * log(i*z + sqrt(1 - z^2))
c_asin = - unit_i * log(unit_i * z + sqrt(1.0d0 - z**2))
end function casin_f
end module complex_acos_asin
program acos_asin_computation
use complex_acos_asin
!use iso_fortran_env
implicit none
complex*8 :: z1, z2, z1_acos, z2_acos, z1_asin, z2_asin
z1 = cmplx(5e8, 5e8, 2)
z2 = cmplx(5e8, -5e8, 2)
z1_acos = cacos_f(z1)
z2_acos = cacos_f(z2)
print *, "cacos(z1) = ", z1_acos
print *, "cacos(z2) = ", z2_acos
print *, ' asin(z1) from acos(z1) = ', 3.14159d0/2-z1_acos
print *, 'acos(z) z=6.7e7-6.7e7i', acos(cmplx(6.7e7, -6.7e7, 2))
z1_asin = casin_f(z1)
z2_asin = casin_f(z2)
print *, "casin(z1) = ", z1_asin
print *, "casin(z2) = ", z2_asin
print *, ' asin(z2) from acos(z2) = ', 3.14159d0/2-z2_acos
end program acos_asin_computation
|
I tested the values in Octave and Scilab and both give the expected solutions
| Code: | z1=5e8 + %i*5e8
acos(z1) = 0.7853982 - 21.069839i
asin(z1) = 0.7853982 + 21.069839i
z1=5e8 - %i*5e8
acos(z2) = 0.7853982 + 21.069839i
asin(z2) = 0.7853982 - 21.069839i |
Is it the sign of the imaginary part causing issues?
Lester |
|
| Back to top |
|
|
Kenneth_Smith
Joined: 18 May 2012
Posts: 797
Location: Hamilton, Lanarkshire, Scotland.
|
| Posted: Sun Mar 02, 2025 2:33 pm Post subject: |
|
|
Thanks Lester,
I think this is what is required:
| Code: | function arcsinCD(z) result (returnval)
integer, parameter :: wp = kind(1.d0)
complex(kind=wp), intent(in) :: z
complex(kind=wp) :: returnval
complex(kind=wp), parameter :: j = cmplx(0.0_wp,1.0_wp,kind=wp)
if (aimag(z) .lt. 0.0_wp) then
returnval = -j*log(sqrt(1.0_wp - z*z) + j*z)
else
returnval = j*log(sqrt(1.0_wp - z*z) - j*z)
end if
end function arcsinCD
function arccosCD(z) result (returnval)
integer, parameter :: wp = kind(1.d0)
complex(kind=wp), intent(in) :: z
complex(kind=wp) :: returnval
complex(kind=wp), parameter :: j = cmplx(0.0_wp,1.0_wp,kind=wp)
complex(kind=wp), parameter :: pi2 = 2.0_wp*atan(1.0_wp)
if (aimag(z) .lt. 0.0_wp) then
returnval = pi2 - arcsinCD(z)
else
returnval = -j*log(j*sqrt( 1.0_wp - z*z) + z)
end if
end function arccosCD |
Here is a test program using these functions:
https://www.dropbox.com/scl/fi/guletuqk37wnocm0mjyk9/test.f95?rlkey=v9ulsz07iqs7enw77gfdy7nsa&st=jwsh8on9&dl=0
And here is the output produced using iFort
https://www.dropbox.com/scl/fi/ilncipneh6d9uh9qjrvr0/log.txt?rlkey=fvktbzfovpvtf5hbqrsn6snvz&st=3ulijujr&dl=0 |
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 146
Location: United Kingdom
|
| Posted: Sun Mar 02, 2025 5:25 pm Post subject: |
|
|
Hi Ken
Good to see you found a solution to the problem.
Lester |
|
| Back to top |
|
|
Kenneth_Smith
Joined: 18 May 2012
Posts: 797
Location: Hamilton, Lanarkshire, Scotland.
|
| Posted: Mon Mar 03, 2025 2:39 pm Post subject: |
|
|
Unfortunately the two functions I posted previously fail to catch the potential underflow to the log function for the case of Z with large (positive or negative) real part and zero imaginary part.
Here are the updated functions:
| Code: | function arcsinCD(z) result (returnval)
integer, parameter :: wp = kind(1.d0)
complex(kind=wp), intent(in) :: z
complex(kind=wp) :: returnval
complex(kind=wp), parameter :: j = cmplx(0.0_wp,1.0_wp,kind=wp)
complex(kind=wp) :: z2, sqrt_term, arglog1, arglog2
if (aimag(z) .lt. 0.0_wp) then
! Imaginary part is less than zero, use eqn 1
returnval = -j*log(sqrt(1.0_wp - z*z) + j*z)
return
else if (aimag(z) .gt. 0.0_wp) then
! Imaginary part is greater than zero, use eqn 2
returnval = j*log(sqrt(1.0_wp - z*z) - j*z)
return
else
! Imaginary part is equal to zero
if ( abs(real(z,kind=wp)) .le. 1.0_wp) then
! If real part is between -1 and 1, the result is not complex
! return the complex number = (asin(real_part), 0.0)
returnval = cmplx(asin(real(z,kind=wp)),0.0_wp,kind=wp)
return
else
! Deal with potential underflow for argument to log function
! for large real part
z2 = z*z
sqrt_term = sqrt(1.0_wp - z2)
arglog1 = sqrt_term + j*z ! depending on sign of the large real
arglog2 = sqrt_term - j*z ! part one of these may underflow to
! (0.0,0.0)
! Use the term with with largest magnitude to avoid (0.0,0.0)
if ( abs (arglog1) .gt. abs(arglog2) ) then
returnval = -j*log(arglog1)
else
returnval = j*log(arglog2)
end if
! Flip sign of imaginary part if negative
if ( aimag(returnval) .lt. 0.0_wp) then
returnval = cmplx(real(returnval,kind=wp), - aimag(returnval), kind=wp)
end if
end if
end if
end function arcsinCD
function arccosCD(z) result (returnval)
integer, parameter :: wp = kind(1.d0)
complex(kind=wp), intent(in) :: z
complex(kind=wp) :: returnval
complex(kind=wp), parameter :: pi2 = 2.0_wp*atan(1.0_wp)
returnval = pi2 - arcsinCD(z)
end function arccosCD
|
Here is my updated test program:
https://www.dropbox.com/scl/fi/guletuqk37wnocm0mjyk9/test.f95?rlkey=v9ulsz07iqs7enw77gfdy7nsa&st=3skgca30&dl=0
And the results obtained using iFort:
https://www.dropbox.com/scl/fi/ilncipneh6d9uh9qjrvr0/log.txt?rlkey=fvktbzfovpvtf5hbqrsn6snvz&st=419qv461&dl=0
No discrepancies between these functions and the native iFort acos(z).
Turned out to be a lot more "complex" than I suspected. |
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 146
Location: United Kingdom
|
| Posted: Wed Mar 05, 2025 2:13 pm Post subject: |
|
|
Hi Ken
I had a look at tweaking the module to cover all domains and it seems to work (so far) with your modifications.
| Code: | module complex_acos_asin
implicit none
interface arcsin
module procedure arcsinR, arcsinCD
end interface arcsin
interface arccos
module procedure arccosR, arccosCD
end interface arccos
contains
function arcsinCD(z) result (returnval)
integer, parameter :: wp = kind(1.d0)
complex(kind=wp), intent(in) :: z
complex(kind=wp) :: returnval
complex(kind=wp), parameter :: j = cmplx(0.0_wp,1.0_wp,kind=wp)
if (aimag(z) < 0.0_wp) then
returnval = -j*log(sqrt(1.0_wp - z*z) + j*z)
else
returnval = j*log(sqrt(1.0_wp - z*z) - j*z)
end if
end function arcsinCD
function arccosCD(z) result (returnval)
integer, parameter :: wp = kind(1.d0)
complex(kind=wp), intent(in) :: z
complex(kind=wp) :: returnval
complex(kind=wp), parameter :: j = cmplx(0.0_wp,1.0_wp,kind=wp)
complex(kind=wp), parameter :: pi2 = 2.0_wp*atan(1.0_wp)
if (aimag(z) < 0.0_wp) then
returnval = pi2 - arcsinCD(z)
else
returnval = -j*log(j*sqrt( 1.0_wp - z*z) + z)
end if
end function arccosCD
function arcsinR(x) result (returnval)
integer, parameter :: wp = kind(1.d0)
real(kind=wp), intent(in) :: x
complex(kind=wp) :: returnval
if (abs(x) <= 1.0_wp) then
returnval = asin(x)
else
returnval = arcsinCD(cmplx(x, 0.0_wp, kind=wp))
end if
end function arcsinR
function arccosR(x) result (returnval)
integer, parameter :: wp = kind(1.d0)
real(kind=wp), intent(in) :: x
complex(kind=wp) :: returnval
if (abs(x) <= 1.0_wp) then
returnval = acos(x)
else
returnval = arccosCD(cmplx(x, 0.0_wp, kind=wp))
end if
end function arccosR
end module complex_acos_asin |
|
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 146
Location: United Kingdom
|
| Posted: Wed Mar 05, 2025 2:21 pm Post subject: |
|
|
Test program to check the results of various inputs
| Code: | program test_asin_acos
use complex_acos_asin
use iso_fortran_env
implicit none
integer :: i
integer, parameter :: wp=kind(1.d0)
real(wp) :: real_part_sin, imag_part_sin, real_part_cos, imag_part_cos
complex(wp) :: asin_res, acos_res, z(2)
! Check that the results of the domain limits are correct for example x = -10 to 10
write( *, '(a)') ' n arcsin(n) arccos(n) '
write( *, '(a)') '=========================================================================='
do i = -10, 10, 1
acos_res = arccos(real(i,kind=wp))
asin_res = arcsin(real(i,kind=wp))
real_part_sin = real(asin_res)
imag_part_sin = aimag(asin_res)
real_part_cos = real(acos_res)
imag_part_cos = aimag(acos_res)
! Output formatted results in normal mathematical style
if (imag_part_sin >= 0.0) then
if (imag_part_cos >= 0.0) then
write(*, '(i5,2x,f15.12, " + ", f12.10, "i", 2x, f15.12, " + ", f13.10, "i")') &
i,real_part_sin, imag_part_sin, real_part_cos, imag_part_cos
else
write(*, '(i5,2x,f15.12, " + ", f12.10, "i", 2x, f15.12, " - ", f13.10, "i")') &
i,real_part_sin, imag_part_sin, real_part_cos, abs(imag_part_cos)
end if
else
if (imag_part_cos >= 0.0) then
write(*, '(i5,2x,f15.12, " - ", f12.10, "i", 2x, f15.12, " + ", f13.10, "i")') &
i,real_part_sin, abs(imag_part_sin), real_part_cos, imag_part_cos
else
write(*, '(i5,2x,f15.12, " - ", f12.10, "i", 2x, f15.12, " - ", f13.10, "i")') &
i,real_part_sin, abs(imag_part_sin), real_part_cos, abs(imag_part_cos)
end if
end if
end do
! verify complex input
z(1) = cmplx(5e8, 5e8, kind=wp)
z(2) = cmplx(5e8, -5e8, kind=wp)
write(*,'()')
print *, 'z1 z2', z
print *, 'z1 (arcsin, arccos = ', arcsin(z(1)), arccos(z(1))
print *, 'z2 (arcsin, arccos = ', arcsin(z(2)), arccos(z(2))
print *, 'arcsin(1.0), arcsin(1.5) ', real(arcsin(1.0_wp)), arcsin(1.5_wp)
print *, arccos(1.0_wp), arccos(0.5_wp)
end program test_asin_acos |
|
|
| Back to top |
|
|
Kenneth_Smith
Joined: 18 May 2012
Posts: 797
Location: Hamilton, Lanarkshire, Scotland.
|
| Posted: Wed Mar 05, 2025 5:54 pm Post subject: |
|
|
Thanks Lester,
This program, which uses the module you posted today, shows the failure which caused me to post the two revised functions on Monday.
| Code: | program test_asin_acos
use complex_acos_asin
use iso_fortran_env
implicit none
integer :: i
integer, parameter :: wp=kind(1.d0)
complex(wp) :: z(2)
! verify complex input
z(1) = cmplx(5e8, 0.0, kind=wp)
z(2) = cmplx(-5e8, 0.0, kind=wp)
write(*,'()')
print *, 'z1 z2', z
print *, 'z1 (arcsin, arccos = ', arcsin(z(1)), arccos(z(1)) ! Fails Re and Im parts to complex log both zero
print *, 'z2 (arcsin, arccos = ', arcsin(z(2)), arccos(z(2)) ! OK
end program test_asin_acos |
Large real part and zero imaginary part causes the argument to complex LOG to be zero.
Here is a MathCad worksheet which shows why this happens:
https://www.dropbox.com/scl/fi/y7am7f8tk48a6zvjteyq9/Screenshot-2025-03-05-164523.jpg?rlkey=wmrsnyqts4yfecs9egyh1fprg&st=td6ims5k&dl=0 |
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 146
Location: United Kingdom
|
| Posted: Wed Mar 05, 2025 6:30 pm Post subject: |
|
|
Thanks for the update Ken. Certainly is an interesting problem.
Lester |
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 146
Location: United Kingdom
|
| Posted: Thu Mar 06, 2025 12:47 pm Post subject: |
|
|
Hi Ken
I updated the code with your latest revisions and it worked fine; the results match the output from Python using cmath as a check.
Interesting point, but Scilab and Octave both give the opposite sign for the imaginary part when testing arcsin(1.5), arccos(1.5).
Thanks again for the reference on this "complex" problem. It is detailed in the NIST Handbook of Mathematical Functions.
Lester |
|
| Back to top |
|
|
Kenneth_Smith
Joined: 18 May 2012
Posts: 797
Location: Hamilton, Lanarkshire, Scotland.
|
| Posted: Fri Mar 07, 2025 12:07 pm Post subject: |
|
|
Lester,
Thanks for the reference. That�s not a text I am familiar with, as there is a copy of Abramowitz and Stegun on my bookshelf � which is just a little older than I am! I�ll track down a copy of this alternative text.
The convention for computing the inverse cosine and inverse sine of values outside the range [−1,1] depends on how the principal branch of the complex logarithm function has been defined. There are two possibilities (both mathematically valid) which is why you see a difference between Scilab/Octave and Python.
At the moment the functions I posted on Wednesday succeed for any complex number on the complex plane within a circle of radius 1.0E+154. Beyond this the z*z term in sqrt(1 � z*z) overflows, since huge(1.0_wp) is approx. 1.0E308, so sqrt(huge(1.0_wp)) is 1.0E+154.
I will look at the possibility of normalising/scaling the sqrt(1 � z*z) calculation to increase the radius of the circle on the complex plane in which the functions succeed. Probably unnecessary, but I�m curious about this. |
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 146
Location: United Kingdom
|
| Posted: Fri Mar 07, 2025 4:25 pm Post subject: |
|
|
Hi Ken
Here is a link to the NIST online mathematical library; it is updated on a regular basis from what I can tell and is a companion to the book (no longer in print).
https://dlmf.nist.gov/
Lester |
|
| Back to top |
|
|
|
|
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum
|
Powered by phpBB © 2001, 2005 phpBB Group
|
FAQ
Search
Memberlist
Usergroups
Register
Profile
Log in to check your private messages
Log in
