Mercurial > mxe-octave

view src/mingw-w64-complex-inverse-trig.patch @ 6356:6f3c099c0d38

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
mingw-w64: Use approximations for casinh if real part is small (bug #62332). * src/mingw-w64-complex-inverse-trig.patch: Use approximations for casinh if real part of input is small and the absolute imaginary part is smaller than 1.
author Markus Mützel <markus.muetzel@gmx.de>
date Sun, 07 Aug 2022 20:09:23 +0200
parents 9d1395334e68
children dc5ad8056086
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

From 98d5b204fe689ac35d348a1742723ec561a2c5f7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Markus=20M=C3=BCtzel?= <markus.muetzel@gmx.de>
Date: Sun, 7 Aug 2022 13:21:29 +0200
Subject: [PATCH 1/2] casinh: Use symmetries also for default case.

---
 mingw-w64-crt/complex/casinh.def.h | 46 ++++++++++++++++--------------
 1 file changed, 25 insertions(+), 21 deletions(-)

diff --git a/mingw-w64-crt/complex/casinh.def.h b/mingw-w64-crt/complex/casinh.def.h
index cce2d864e..d309e6392 100644
--- a/mingw-w64-crt/complex/casinh.def.h
+++ b/mingw-w64-crt/complex/casinh.def.h
@@ -89,32 +89,36 @@ __FLT_ABI(casinh) (__FLT_TYPE __complex__ 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)
+  /* Use symmetries to perform the calculation in the first quadrant. */
+  __FLT_TYPE arz = __FLT_ABI(fabs) (__real__ z);
+  __FLT_TYPE aiz = __FLT_ABI(fabs) (__imag__ z);
+
+  if (arz >= __FLT_CST(1.0)/__FLT_EPSILON
+      || aiz >= __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;
+    Use that approximation to avoid overflow when squaring. */
+    __real__ x = arz;
+    __imag__ x = aiz;
+    ret = __FLT_ABI(clog) (x);
+    __real__ ret += M_LN2;
   }
+  else
+  {
+    __real__ x = (arz - aiz) * (arz + aiz) + __FLT_CST(1.0);
+    __imag__ x = __FLT_CST(2.0) * arz * aiz;
+
+    x = __FLT_ABI(csqrt) (x);
 
-  __real__ x = (__real__ z - __imag__ z) * (__real__ z + __imag__ z) + __FLT_CST(1.0);
-  __imag__ x = __FLT_CST(2.0) * __real__ z * __imag__ z;
+    __real__ x += arz;
+    __imag__ x += aiz;
 
-  x = __FLT_ABI(csqrt) (x);
+    ret = __FLT_ABI(clog) (x);
+  }
 
-  __real__ x += __real__ z;
-  __imag__ x += __imag__ z;
+  /* adjust signs for input quadrant */
+  __real__ ret = __FLT_ABI(copysign) (__real__ ret, __real__ z);
+  __imag__ ret = __FLT_ABI(copysign) (__imag__ ret, __imag__ z);
 
-  return __FLT_ABI(clog) (x);
+  return ret;
 }
-- 
2.35.3.windows.1

From 099f1cd2d4261573d584aa415f47742fa977c773 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Markus=20M=C3=BCtzel?= <markus.muetzel@gmx.de>
Date: Sun, 7 Aug 2022 13:23:02 +0200
Subject: [PATCH 2/2] casinh: Use approximation for small real part and
 absolute imaginary part < 1

---
 mingw-w64-crt/complex/casinh.def.h | 31 ++++++++++++++++++++++++++++++
 1 file changed, 31 insertions(+)

From 575a070e3fb360fced10a458e72a298d3071c42e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Markus=20M=C3=BCtzel?= <markus.muetzel@gmx.de>
Date: Sun, 7 Aug 2022 19:05:21 +0200
Subject: [PATCH 2/2] casinh: Use approximations for small real part and
 absolute imaginary part < 1

---
 mingw-w64-crt/complex/casinh.def.h | 31 ++++++++++++++++++++++++++++++
 1 file changed, 31 insertions(+)

diff --git a/mingw-w64-crt/complex/casinh.def.h b/mingw-w64-crt/complex/casinh.def.h
index d309e6392..a5447bfcd 100644
--- a/mingw-w64-crt/complex/casinh.def.h
+++ b/mingw-w64-crt/complex/casinh.def.h
@@ -103,6 +103,37 @@ __FLT_ABI(casinh) (__FLT_TYPE __complex__ z)
     ret = __FLT_ABI(clog) (x);
     __real__ ret += M_LN2;
   }
+  else if (arz*arz <= __FLT_EPSILON && aiz < __FLT_CST(1.0))
+  {
+    /* Taylor series expansion around arz=0 for z + sqrt(z*z + 1):
+    c = arz + sqrt(1-aiz^2) + arz^2 / (2*(1-aiz^2)^(3/2)) + i*(aiz + arz*aiz / sqrt(1-aiz^2)) + O(arz^4)
+    Identity: clog(c) = log(|c|) + i*arg(c)
+    For real part of result:
+    |c| = 1 + arz / sqrt(1-aiz^2) + arz^2/(2*(1-aiz^2)) + O(arz^3)  (Taylor series expansion)
+    For imaginary part of result:
+    c = 1/sqrt(1-aiz^2) * ((1-aiz^2) + arz*sqrt(1-aiz^2) + arz^2/(2*(1-aiz^2)) + i*aiz*(sqrt(1-aiz^2)+arz)) + O(arz^3)
+    */
+    __FLT_TYPE onemaiz2 = (__FLT_CST(1.0)+aiz)*(__FLT_CST(1.0)-aiz);
+    __FLT_TYPE s1maiz2 = __FLT_ABI(sqrt) (onemaiz2);
+    __FLT_TYPE arz2red = arz * arz / __FLT_CST(2.0) / s1maiz2;
+    __real__ ret = __FLT_ABI(log1p) ((arz + arz2red) / s1maiz2);
+    __imag__ ret = __FLT_ABI(atan2) (aiz * (s1maiz2 + arz),
+                                     onemaiz2 + arz*s1maiz2 + arz2red);
+  }
+  else if (arz <= __FLT_EPSILON && aiz < __FLT_CST(1.0))
+  {
+    /* Taylor series expansion around arz=0 for z + sqrt(z*z + 1):
+    c = arz + sqrt(1-aiz^2) + i*(aiz + arz*aiz / sqrt(1-aiz^2)) + O(arz^2)
+    Identity: clog(c) = log(|c|) + i*arg(c)
+    For real part of result:
+    |c| = 1 + arz / sqrt(1-aiz^2) + O(arz^2)  (Taylor series expansion)
+    For imaginary part of result:
+    c = (arz + sqrt(1-aiz^2))/sqrt(1-aiz^2) * (sqrt(1-aiz^2) + i*aiz) + O(arz^6)
+    */
+    __FLT_TYPE s1maiz2 = __FLT_ABI(sqrt) ((__FLT_CST(1.0)+aiz)*(__FLT_CST(1.0)-aiz));
+    __real__ ret = __FLT_ABI(log1p) (arz / s1maiz2);
+    __imag__ ret = __FLT_ABI(atan2) (aiz, s1maiz2);
+  }
   else
   {
     __real__ x = (arz - aiz) * (arz + aiz) + __FLT_CST(1.0);
-- 
2.35.3.windows.1