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

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

Add a New Comment (Rich Markup)
   

Comment Type & Canned Response:
       

 

( Jump to the original submission )

Mon 29 Aug 2022 10:46:39 AM UTC, comment #11: 

The last patch has been accepted upstream:
https://sourceforge.net/p/mingw-w64/mingw-w64/ci/dbeb7b2720824e10292c9356fe7b0da89c6b3b2c/

Closing as fixed.

Markus Mützel <mmuetzel>
Project Member
Sun 28 Aug 2022 03:10:15 PM UTC, comment #10: 
Markus Mützel <mmuetzel>
Project Member
Sun 14 Aug 2022 11:35:11 AM UTC, comment #9: 
Markus Mützel <mmuetzel>
Project Member
Mon 08 Aug 2022 06:28:27 PM UTC, comment #8: 

I updated the patch in MXE Octave to improve the accuracy of catanh for small real part (and to fix a thinko in the patch for casinh):
https://hg.octave.org/mxe-octave/rev/dc5ad8056086

With it, I see the following for the example from comment #0 on Windows:

>> for ij=1:20, x = complex(1/sqrt(2), 10^(-3*ij)); printf('%.16e%+.16ei -> %.16e%+.16ei\n', real(x), imag(x), real(acos(x)),imag(acos(x))); end
7.0710678118654746e-01+1.0000000000000002e-03i -> 7.8539916339344829e-01-1.4142116767647246e-03i
7.0710678118654746e-01+1.0000000000000006e-06i -> 7.8539816339844837e-01-1.4142135623447793e-06i
7.0710678118654746e-01+1.0000000000000009e-09i -> 7.8539816339744850e-01-1.4142135623730963e-09i
7.0710678118654746e-01+1.0000000000000014e-12i -> 7.8539816339744839e-01-1.4142135623730968e-12i
7.0710678118654746e-01+1.0000000000000017e-15i -> 7.8539816339744839e-01-1.4142135623730974e-15i
7.0710678118654746e-01+1.0000000000000018e-18i -> 7.8539816339744828e-01-1.4142135623730976e-18i
7.0710678118654746e-01+1.0000000000000024e-21i -> 7.8539816339744828e-01-1.4142135623730982e-21i
7.0710678118654746e-01+1.0000000000000025e-24i -> 7.8539816339744828e-01-1.4142135623730985e-24i
7.0710678118654746e-01+1.0000000000000027e-27i -> 7.8539816339744828e-01-1.4142135623730988e-27i
7.0710678118654746e-01+1.0000000000000031e-30i -> 7.8539816339744828e-01-1.4142135623730993e-30i
7.0710678118654746e-01+1.0000000000000035e-33i -> 7.8539816339744828e-01-1.4142135623730999e-33i
7.0710678118654746e-01+1.0000000000000040e-36i -> 7.8539816339744828e-01-1.4142135623731006e-36i
7.0710678118654746e-01+1.0000000000000042e-39i -> 7.8539816339744828e-01-1.4142135623731009e-39i
7.0710678118654746e-01+1.0000000000000045e-42i -> 7.8539816339744828e-01-1.4142135623731013e-42i
7.0710678118654746e-01+1.0000000000000050e-45i -> 7.8539816339744828e-01-1.4142135623731019e-45i
7.0710678118654746e-01+1.0000000000000051e-48i -> 7.8539816339744828e-01-1.4142135623731023e-48i
7.0710678118654746e-01+1.0000000000000053e-51i -> 7.8539816339744828e-01-1.4142135623731026e-51i
7.0710678118654746e-01+1.0000000000000058e-54i -> 7.8539816339744828e-01-1.4142135623731032e-54i
7.0710678118654746e-01+1.0000000000000060e-57i -> 7.8539816339744828e-01-1.4142135623731035e-57i
7.0710678118654746e-01+1.0000000000000066e-60i -> 7.8539816339744828e-01-1.4142135623731042e-60i

Marking as ready for test.

Markus Mützel <mmuetzel>
Project Member
Sun 07 Aug 2022 06:24:50 PM UTC, comment #7: 

Octave uses the math functions from the standard libraries.
It looks like the complex inverse trigonometric and hyperbolic functions of mingw-w64 have accuracy issues for some regimes.
`cacos` and `casin` both depend on the implementation of `casinh`.
I pushed a patch to MXE Octave that should improve the accuracy of these function if the absolute value of the imaginary part of the input to `casinh` is smaller than 1 and the absolute value of the real part is small:
https://hg.octave.org/mxe-octave/rev/6f3c099c0d38

That should hopefully cover the case from comment #0.

I didn't have a look at the implementation of `catan` (respectively, `catanh`) yet.

Markus Mützel <mmuetzel>
Project Member
Wed 06 Jul 2022 06:08:09 AM UTC, comment #6: 

Typo: "internet" should read "integer". (Curse you, touch keyboard 😜)

Markus Mützel <mmuetzel>
Project Member
Wed 06 Jul 2022 06:06:43 AM UTC, comment #5: 

I rebuilt Octave with the updated MinGW CRT (in MSYS2). With that, the precision of powers (with negative internet exponent) is indeed better than it was before. However, it doesn't make a difference for the test cases here.

Markus Mützel <mmuetzel>
Project Member
Mon 04 Jul 2022 02:03:30 PM UTC, comment #4: 

A change to the implementation of the power functions landed in upstream MinGW:
https://sourceforge.net/p/mingw-w64/mailman/message/37668190/

I'll try and check if that makes a difference here when I come around to it. (And possible cherry pick for MXE Octave?)

Markus Mützel <mmuetzel>
Project Member
Mon 27 Jun 2022 06:47:15 AM UTC, comment #3: 

I'm not sure if this is related. But the implementation of the power functions in MinGW seems to have accuracy issues for some input. Maybe, the functions here are implemented as some kind of power series. If that is the case, fixing the accuracy issues upstream might also help here.

See also this thread on the mingw-w64 mailing list:
https://sourceforge.net/p/mingw-w64/mailman/mingw-w64-public/thread/2ad75ebab0eb95e7612936debc11445c%40m-labs.hk/#msg37663906

Markus Mützel <mmuetzel>
Project Member
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 jwe (Updated the item)
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -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.

     

     

     

     

    Follow 7 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-08-29 mmuetzel StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2022-08-08 mmuetzel StatusIn Progress Ready For Test
    2022-08-07 mmuetzel StatusNone In Progress
        Release7.1.0 dev
    2022-08-03 jwe CategoryNone Libraries
    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