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. Markus
From de5c54517f2102dc43d54745b307e8c00a27000e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Markus=20M=C3=BCtzel?= <[email protected]> Date: Tue, 2 Aug 2022 17:07:11 +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..1f03befae 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 plane. */ + __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
