Mercurial > mxe-octave
view src/mingw-w64-complex-inverse-trig.patch @ 6354:9d1395334e68
mingw-w64: Handle negative zero input for cacosh (bug #49091).
* src/mingw-w64-complex-inverse-trig.patch: Handle negative zero input for
cacosh.
author | Markus Mützel <markus.muetzel@gmx.de> |
---|---|
date | Fri, 05 Aug 2022 19:32:03 +0200 |
parents | 0825abaf61a7 |
children | 6f3c099c0d38 |
line wrap: on
line source
From 59fd6e6ab7c93545175fa6ff46532b7a80dbdbfb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Markus=20M=C3=BCtzel?= <markus.muetzel@gmx.de> 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?= <markus.muetzel@gmx.de> 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 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 From 9fcac91fc50120a5e6acfa94ba0a917c7faf2306 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Markus=20M=C3=BCtzel?= <markus.muetzel@gmx.de> Date: Fri, 5 Aug 2022 18:36:49 +0200 Subject: [PATCH] cacosh: Handle negative zero input --- mingw-w64-crt/complex/cacosh.def.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mingw-w64-crt/complex/cacosh.def.h b/mingw-w64-crt/complex/cacosh.def.h index 1f03befae..addf2958e 100644 --- a/mingw-w64-crt/complex/cacosh.def.h +++ b/mingw-w64-crt/complex/cacosh.def.h @@ -106,7 +106,7 @@ __FLT_ABI(cacosh) (__FLT_TYPE __complex__ z) x = __FLT_ABI(csqrt) (x); - if (__real__ z < __FLT_CST(0.0)) + if (__FLT_ABI(signbit) (__real__ z)) x = -x; __real__ x += __real__ z; @@ -114,7 +114,7 @@ __FLT_ABI(cacosh) (__FLT_TYPE __complex__ z) ret = __FLT_ABI(clog) (x); - if (__real__ ret < __FLT_CST(0.0)) + if (__FLT_ABI(signbit) (__real__ ret)) ret = -ret; return ret; -- 2.35.3.windows.1