changeset 6335:0825abaf61a7

mingw-w64: Special case for large input for complex inverse trigonometric functions (bug #49091). * src/mingw-w64-complex-inverse-trig.patch: Add patch for mingw-w64 CRT that implements a special case for large input for complex inverse trigonometric functions cacosh and casinh (and therefore functions that are derived from it, like cacos and casin). * dist-files.mk: Add new patch to list.
author Markus Mützel <markus.muetzel@gmx.de>
date Sat, 23 Jul 2022 19:00:09 +0200
parents 892ce68bc32f
children 950580a30442
files dist-files.mk src/mingw-w64-complex-inverse-trig.patch
diffstat 2 files changed, 89 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- a/dist-files.mk	Fri Jul 22 15:05:59 2022 -0400
+++ b/dist-files.mk	Sat Jul 23 19:00:09 2022 +0200
@@ -352,6 +352,7 @@
   mingw-texinfo-1-fixes.patch \
   mingw-texinfo-2-makeinfo-non-ASCII-perl.patch \
   mingw-utils-1-portability-fix.patch \
+  mingw-w64-complex-inverse-trig.patch \
   mingw-w64.mk \
   mingw-zeromq-1-fixes.patch \
   mingwrt.mk \
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/mingw-w64-complex-inverse-trig.patch	Sat Jul 23 19:00:09 2022 +0200
@@ -0,0 +1,88 @@
+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 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
+