Puzzle asin(z) and acos(z)

The complex arcsin and arccos functions can each be expressed in two ways:

    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 – zz is common to three of these, and the addition of +jz 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:

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:

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).

Continued .....

For the full code above, with the lines problematic to FTN not commented out, we get:

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.

Ken

I have made a note to look into this.

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)

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

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!

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:

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

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

Thanks Lester,

I think this is what is required:

  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

Hi Ken

Good to see you found a solution to the problem.

Lester

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:

  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.

Hi Ken

I had a look at tweaking the module to cover all domains and it seems to work (so far) with your modifications.

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

Test program to check the results of various inputs

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

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.

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

Thanks for the update Ken. Certainly is an interesting problem.

Lester

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

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 zz term in sqrt(1 – zz) 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.

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

Many thanks for this valuable analysis. I have implemented Ken's method in the Win32/x64 FTN95 library for the next release.