bugGNU Octave - Bugs: bug #62332, [MinGW] acos(z), asin(z) and...

 
 

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

Submitted by:  None
Submitted on:  Mon 18 Apr 2022 09:41:59 PM UTC
 
Category:  None Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  None Assigned to:  None
Originator Name:  Cristiano Originator Email:  -email is unavailable-
Open/Closed:  * Open Release:  * 7.1.0
Operating System:  * Microsoft Windows
* Mandatory Fields

Add a New Comment (Rich Markup)
   

Comment Type & Canned Response:
       

 

Mon 18 Apr 2022 11:35:21 PM UTC, comment #2: 

Yes Rik, this summarizes the problem perfectly.

comentário #1:

> Is the summary that the C++ math library on MinGW is inaccurate?
>
> It seems that a test for the behavior can be boiled down to
>

> imag (acos (complex (cosd (45), 1e-18)))

>
> This should return -1.4e-18, but returns 0 on MinGW.  Is that correct?

Anonymous
Mon 18 Apr 2022 10:43:29 PM UTC, comment #1: 

Is the summary that the C++ math library on MinGW is inaccurate?

It seems that a test for the behavior can be boiled down to

imag (acos (complex (cosd (45), 1e-18)))

This should return -1.4e-18, but returns 0 on MinGW.  Is that correct?

Rik <rik5>
Project Administrator
Mon 18 Apr 2022 09:41:59 PM UTC, original submission:  

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.

Anonymous

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by None (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only project members can vote.

     

     

     

     

    Follows 1 latest change.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-04-18 rik5 Summaryacos(z), asin(z) and atan(z) , ( z = x + yi ) return wrong result for imaginary part lower than 1e-12 => [MinGW] acos(z), asin(z) and atan(z) , ( z = x + yi ) return wrong result for imaginary part lower than 1e-12

    Back to the top


    Powered by Savane 3.9