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

Group

Main
- Main
- View members
- Search
Homepage
Download
Mailing lists
Source code
Bugs
- Submit new
- Browse
- Reset to open
- Digest
- Dependencies
- ~~Export~~
- Get statistics
- Search
Tasks
- ~~Submit new~~
- Browse
- Reset to open
- Digest
- Dependencies
- ~~Export~~
- Get statistics
- Search
Patches
- Submit new
- Browse
- Reset to open
- Digest
- Dependencies
- ~~Export~~
- Get statistics
- Search

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

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>
Sun 28 Aug 2022 03:10:15 PM UTC, comment #10: The first couple of patches have been accepted upstream: https://sourceforge.net/p/mingw-w64/mingw-w64/ci/3ba7dce871477ba9739d5f1e00ccec4cc1ccc693/ https://sourceforge.net/p/mingw-w64/mingw-w64/ci/1dd2a4993b4d2431e4ab2320d3558232f9b2ef45/ I sent the last patch (for catanh) for review to the mingw64 mailing list here: https://sourceforge.net/p/mingw-w64/mailman/mingw-w64-public/thread/trinity-5263ce84-db46-481d-a408-80caf3f90fb9-1661695236627%403c-app-gmx-bs26/#msg37699549	Markus Mützel <mmuetzel>
Sun 14 Aug 2022 11:35:11 AM UTC, comment #9: I sent the first couple of patches (for asinh) to the mingw64 mailing list: https://sourceforge.net/p/mingw-w64/mailman/mingw-w64-public/thread/trinity-393c5a07-b7ab-43ba-8cfa-bb6bfec2cbbb-1660476467217%403c-app-gmx-bap70/#msg37693407	Markus Mützel <mmuetzel>
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>
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>
Wed 06 Jul 2022 06:08:09 AM UTC, comment #6: Typo: "internet" should read "integer". (Curse you, touch keyboard 😜)	Markus Mützel <mmuetzel>
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>
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>
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>
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>
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^(-3ij) ); 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^(-3ij) ); 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^(-3ij) ); 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 group members can vote.

Please enter the title of George Orwell's famous dystopian book (it's a date):

Follow 7 latest changes.

Date	Changed by	Updated Field	Previous Value	Replaced by
2022-08-29	mmuetzel	Status	Ready For Test	Fixed
		Open/Closed	Open	Closed
2022-08-08	mmuetzel	Status	In Progress	Ready For Test
2022-08-07	mmuetzel	Status	None	In Progress
		Release	7.1.0	dev
2022-08-03	jwe	Category	None	Libraries
2022-04-18	rik5	Summary	acos(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

GNU 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

<a href="#postcomment">Post a Comment</a>

<a href="#discussion">Discussion</a>

<a href="#attached">Attached Files</a>

<a href="#dependencies">Dependencies</a>

<a href="#cc">Mail Notification Carbon-Copy List</a>

<a href="#votes">Votes</a>

<a href="#history">History</a>