[Octave-bug-tracker] [bug #62332] acos(z), asin(z) and atan(z) , ( z = x + yi ) return wrong result for imaginary part lower than 1e-12

[anônimo]

URL:
  <https://savannah.gnu.org/bugs/?62332>

                 Summary: acos(z), asin(z) and atan(z) , ( z = x + yi ) return
wrong result for imaginary part lower than 1e-12
                 Project: GNU Octave
            Submitted by: None
            Submitted on: Seg 18 Abr 2022 21:41:59 UTC
                Category: None
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: Incorrect Result
                  Status: None
             Assigned to: None
         Originator Name: Cristiano
        Originator Email: cristiano.ubessi@ufrgs.br
             Open/Closed: Open
                 Release: 7.1.0
         Discussion Lock: Any
        Operating System: Microsoft Windows

    _______________________________________________________

Details:

When trying to calculate derivatives with the complex-step method, i.e.,
appending a small imaginary part to x, i.e.: z= x + yi, with y<1e-16, the arc
cosine or cos^-1 function returns a null imaginary part, which was not
expected to occur.


- Apparently, this has to do with the MinGW distribution. (C++ code also
yields the same results).

- Looking around on the web I have only found a similar report in NumPy, in
this link:
https://github.com/numpy/numpy/issues/6081

Test results:

On Octave 7.1.0 on MS Windows 10 (same on 5.2.0 and 6.4.0):


>> disp( version );computer;for ij=1:20;x=complex( cosd(45) , 10^(-3*ij) );
printf('%e%+ei -> %e%+ei \n', real(x),imag(x),real(acos( x
)),imag(acos(x)) ) ;end
7.1.0
x86_64-w64-mingw32
7.071068e-01+1.000000e-03i -> 7.853992e-01-1.414212e-03i
7.071068e-01+1.000000e-06i -> 7.853982e-01-1.414214e-06i
7.071068e-01+1.000000e-09i -> 7.853982e-01-1.414214e-09i
7.071068e-01+1.000000e-12i -> 7.853982e-01-1.414202e-12i
7.071068e-01+1.000000e-15i -> 7.853982e-01-1.443290e-15i
7.071068e-01+1.000000e-18i -> 7.853982e-01+0.000000e+00i
7.071068e-01+1.000000e-21i -> 7.853982e-01+0.000000e+00i
7.071068e-01+1.000000e-24i -> 7.853982e-01+0.000000e+00i
7.071068e-01+1.000000e-27i -> 7.853982e-01+0.000000e+00i
7.071068e-01+1.000000e-30i -> 7.853982e-01+0.000000e+00i
7.071068e-01+1.000000e-33i -> 7.853982e-01+0.000000e+00i
7.071068e-01+1.000000e-36i -> 7.853982e-01+0.000000e+00i
7.071068e-01+1.000000e-39i -> 7.853982e-01+0.000000e+00i
7.071068e-01+1.000000e-42i -> 7.853982e-01+0.000000e+00i
7.071068e-01+1.000000e-45i -> 7.853982e-01+0.000000e+00i
7.071068e-01+1.000000e-48i -> 7.853982e-01+0.000000e+00i
7.071068e-01+1.000000e-51i -> 7.853982e-01+0.000000e+00i
7.071068e-01+1.000000e-54i -> 7.853982e-01+0.000000e+00i
7.071068e-01+1.000000e-57i -> 7.853982e-01+0.000000e+00i
7.071068e-01+1.000000e-60i -> 7.853982e-01+0.000000e+00i




- Tested on Debian, Octave 6.2.0, on wsl2, and the result is correct (same as
Matlab R2012a).


>> disp( version );computer;for ij=1:20;x=complex( cosd(45) , 10^(-3*ij) );
printf('%e%+ei -> %e%+ei \n', real(x),imag(x),real(acos( x
)),imag(acos(x)) ) ;end
6.2.0
x86_64-pc-linux-gnu
7.071068e-01+1.000000e-03i -> 7.853992e-01-1.414212e-03i
7.071068e-01+1.000000e-06i -> 7.853982e-01-1.414214e-06i
7.071068e-01+1.000000e-09i -> 7.853982e-01-1.414214e-09i
7.071068e-01+1.000000e-12i -> 7.853982e-01-1.414214e-12i
7.071068e-01+1.000000e-15i -> 7.853982e-01-1.414214e-15i
7.071068e-01+1.000000e-18i -> 7.853982e-01-1.414214e-18i
7.071068e-01+1.000000e-21i -> 7.853982e-01-1.414214e-21i
7.071068e-01+1.000000e-24i -> 7.853982e-01-1.414214e-24i
7.071068e-01+1.000000e-27i -> 7.853982e-01-1.414214e-27i
7.071068e-01+1.000000e-30i -> 7.853982e-01-1.414214e-30i
7.071068e-01+1.000000e-33i -> 7.853982e-01-1.414214e-33i
7.071068e-01+1.000000e-36i -> 7.853982e-01-1.414214e-36i
7.071068e-01+1.000000e-39i -> 7.853982e-01-1.414214e-39i
7.071068e-01+1.000000e-42i -> 7.853982e-01-1.414214e-42i
7.071068e-01+1.000000e-45i -> 7.853982e-01-1.414214e-45i
7.071068e-01+1.000000e-48i -> 7.853982e-01-1.414214e-48i
7.071068e-01+1.000000e-51i -> 7.853982e-01-1.414214e-51i
7.071068e-01+1.000000e-54i -> 7.853982e-01-1.414214e-54i
7.071068e-01+1.000000e-57i -> 7.853982e-01-1.414214e-57i
7.071068e-01+1.000000e-60i -> 7.853982e-01-1.414214e-60i



The C++ code below shows the same behaviour when compiled via octave command
line or using the 'bash shell' shipped with octave for Windows:


#include <iostream>
#include <complex>
int main()
{
    typedef typename std::complex<double> Complex;
    Complex x = Complex( 0,0 ) ;
    Complex y = Complex(-1.904e-01, +0);
    Complex z = Complex(-1.904e-01 ,- 9.518e-80 )  ;
    Complex rho = sqrt(x * x + y * y + z * z);
    std::cout << "phi ="<< acos( z / rho ) << std::endl ;
}




- This error may affect only a few people that use the complex-step method
(similarly to automatic differentiation) in code that performs transformation
from a cartesian to a spherical coordinate system. These few people may
already have to implement a complex atan2(y,x) function. Otherwise, the cos
and sin functions are working correctly.


A possible palliative solution for .oct files is to replace the std::acos() by
the following function:


typedef typename std::complex<double> Complex;
Complex acos( Complex & X )
{
      return Complex( acos( X.real() )  , 
             // this is the derivative of cos^-1(x), times the imaginary part
of x                        
             - X.imag() / sqrt( 1 - pow( X.real() , 2 ) ) 
             ) ;
}


Below is the obtained result by using the palliative solution, derivative os
cos^-1 for the imaginary part. This is only useful for the complex-step method
as for large imaginary parts the real part is not affected and the result is
actually wrong.


>> disp( version );computer;for ij=1:20;x=complex( cosd(45) , 10^(-3*ij) );
printf('%e%+ei -> %e%+ei \n', real(x),imag(x), acos(real(x) ),
-imag(x)/sqrt(1-real(x)^2) ) ;end
7.1.0
x86_64-w64-mingw32
7.071068e-01+1.000000e-03i -> 7.853982e-01-1.414214e-03i < Wrong real part
7.071068e-01+1.000000e-06i -> 7.853982e-01-1.414214e-06i
7.071068e-01+1.000000e-09i -> 7.853982e-01-1.414214e-09i
7.071068e-01+1.000000e-12i -> 7.853982e-01-1.414214e-12i
7.071068e-01+1.000000e-15i -> 7.853982e-01-1.414214e-15i
7.071068e-01+1.000000e-18i -> 7.853982e-01-1.414214e-18i
7.071068e-01+1.000000e-21i -> 7.853982e-01-1.414214e-21i
7.071068e-01+1.000000e-24i -> 7.853982e-01-1.414214e-24i
7.071068e-01+1.000000e-27i -> 7.853982e-01-1.414214e-27i
7.071068e-01+1.000000e-30i -> 7.853982e-01-1.414214e-30i
7.071068e-01+1.000000e-33i -> 7.853982e-01-1.414214e-33i
7.071068e-01+1.000000e-36i -> 7.853982e-01-1.414214e-36i
7.071068e-01+1.000000e-39i -> 7.853982e-01-1.414214e-39i
7.071068e-01+1.000000e-42i -> 7.853982e-01-1.414214e-42i
7.071068e-01+1.000000e-45i -> 7.853982e-01-1.414214e-45i
7.071068e-01+1.000000e-48i -> 7.853982e-01-1.414214e-48i
7.071068e-01+1.000000e-51i -> 7.853982e-01-1.414214e-51i
7.071068e-01+1.000000e-54i -> 7.853982e-01-1.414214e-54i
7.071068e-01+1.000000e-57i -> 7.853982e-01-1.414214e-57i
7.071068e-01+1.000000e-60i -> 7.853982e-01-1.414214e-60i



Kind regards,
Cristiano.
*This is my first bug report, any suggestions are welcome.





    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?62332>

_______________________________________________
  Mensagem enviada pelo Savannah
  https://savannah.gnu.org/