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. I'm looking forward to any feedback and hope that something similar could be accepted. Markus
From 59fd6e6ab7c93545175fa6ff46532b7a80dbdbfb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Markus=20M=C3=BCtzel?= <[email protected]> Date: Sat, 23 Jul 2022 15:30:46 +0200 Subject: [PATCH 1/2] casinh: Use approximation for large input. --- mingw-w64-crt/complex/casinh.def.h | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/mingw-w64-crt/complex/casinh.def.h b/mingw-w64-crt/complex/casinh.def.h index 050d885a0..cce2d864e 100644 --- a/mingw-w64-crt/complex/casinh.def.h +++ b/mingw-w64-crt/complex/casinh.def.h @@ -87,6 +87,27 @@ __FLT_ABI(casinh) (__FLT_TYPE __complex__ z) if (r_class == FP_ZERO && i_class == FP_ZERO) return z; + /* casinh(z) = log(z + sqrt(z*z + 1)) */ + + if (__FLT_ABI(fabs) (__real__ z) >= __FLT_CST(1.0)/__FLT_EPSILON + || __FLT_ABI(fabs) (__imag__ z) >= __FLT_CST(1.0)/__FLT_EPSILON) + { + /* For large z, z + sqrt(z*z + 1) is approximately 2*z. + Use that approximation to avoid overflow when squaring. + Additionally, use symmetries to perform the calculation in the first + quadrant. */ + __real__ x = __FLT_ABI(fabs) (__real__ z); + __imag__ x = __FLT_ABI(fabs) (__imag__ z); + x = __FLT_ABI(clog) (x); + __real__ x += M_LN2; + + /* adjust signs for input quadrant */ + __real__ ret = __FLT_ABI(copysign) (__real__ x, __real__ z); + __imag__ ret = __FLT_ABI(copysign) (__imag__ x, __imag__ z); + + return ret; + } + __real__ x = (__real__ z - __imag__ z) * (__real__ z + __imag__ z) + __FLT_CST(1.0); __imag__ x = __FLT_CST(2.0) * __real__ z * __imag__ z; -- 2.35.3.windows.1
From 4b5a229573ca98a622cc15d20bd561dbaec9bbfc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Markus=20M=C3=BCtzel?= <[email protected]> Date: Sat, 23 Jul 2022 18:01:29 +0200 Subject: [PATCH 2/2] cacosh: Use approximation for large input. --- mingw-w64-crt/complex/cacosh.def.h | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/mingw-w64-crt/complex/cacosh.def.h b/mingw-w64-crt/complex/cacosh.def.h index f4ea2da07..a3edd5288 100644 --- a/mingw-w64-crt/complex/cacosh.def.h +++ b/mingw-w64-crt/complex/cacosh.def.h @@ -80,6 +80,27 @@ __FLT_ABI(cacosh) (__FLT_TYPE __complex__ z) return ret; } + /* cacosh(z) = log(z + sqrt(z*z - 1)) */ + + if (__FLT_ABI(fabs) (__real__ z) >= __FLT_CST(1.0)/__FLT_EPSILON + || __FLT_ABI(fabs) (__imag__ z) >= __FLT_CST(1.0)/__FLT_EPSILON) + { + /* For large z, z + sqrt(z*z - 1) is approximately 2*z. + Use that approximation to avoid overflow when squaring. + Additionally, use symmetries to perform the calculation in the positive half + plain. */ + __real__ x = __real__ z; + __imag__ x = __FLT_ABI(fabs) (__imag__ z); + x = __FLT_ABI(clog) (x); + __real__ x += M_LN2; + + /* adjust signs for input */ + __real__ ret = __real__ x; + __imag__ ret = __FLT_ABI(copysign) (__imag__ x, __imag__ z); + + return ret; + } + __real__ x = (__real__ z - __imag__ z) * (__real__ z + __imag__ z) - __FLT_CST(1.0); __imag__ x = __FLT_CST(2.0) * __real__ z * __imag__ z; -- 2.35.3.windows.1
_______________________________________________ Mingw-w64-public mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/mingw-w64-public
