On Tue, 2 Aug 2022, Markus Mützel wrote:
Am 02. August 2022 um 15:09 Uhr schrieb "Martin Storsjö":
On Tue, 26 Jul 2022, Markus Mützel wrote:
Hello,
We noticed in downstream Octave that the results of complex inverse
trigonometric functions are off (sometimes by pi/2) for large input on MinGW.
See also this downstream report:
https://savannah.gnu.org/bugs/index.php?49091
Afaict, this is caused by the following:
The implementation of cacos uses cacosh, and casin uses casinh.
"cacosh" is implemented using the identity: cacosh(z) = log(z + sqrt(z*z - 1))
"casinh" is implemented using the identity: casinh(z) = log(z + sqrt(z*z + 1))
IIUC, "z*z" might overflow or become numerically unstable for large floating point
"z" in both cases.
The proposed patches use approximations for large z that avoid squaring "z" for large values of "z". For large
"z", "z + sqrt(z*z + 1)" or "z + sqrt(z*z - 1)" is approximately "2*z". Additionally, it uses
symmetries to get the correct results for input in any quadrant.
I compared the results for some large z (z=±1e150; z=±1e150i; z=±1e150±1e150i;
for cacos, cacosh, casin, and casinh) to the results of Wolfram Alfa.
E.g.: https://www.wolframalpha.com/input?i=acos%281e150%29
With the proposed approximations, they coincide to 16 significant decimal
digits in the real and imaginary parts for double precision floating point.
Without the approximations (the current base line), the deviations are about
pi/2 for some (e.g., the real part of cacos(1e150+0i) or cacos(-1e150+0i)).
This is the first time I'm writing to this list. So, please let me know if this
is not the correct place or format to propose patches.
This is the right place for such patches.
The patches look good to me (although it took some time to work through
the math - it's been a while since I used my complex math); I didn't try
to check how the symmetries that are used work out here, but it seems
reasonable and I presume you tested those aspects.
Thank you for reviewing.
The symmetries should be correct in theory. And I also tested them with a
couple of different inputs.
Minor nitpick; in cacosh in the comment, "half plain" should be "half
plane".
I'm attaching an updated patch with that spelling error corrected.
I can push the patch later if there's no other disagreement about it.
Looking forward to that.
Pushed now, thanks for the patch!
// Martin
_______________________________________________
Mingw-w64-public mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/mingw-w64-public
