changeset 6357:dc5ad8056086

mingw-w64: Use approximations for catanh if real part is small (bug #62332). * src/mingw-w64-complex-inverse-trig.patch: Use approximations for catanh if real part of input is small. Fix thinko in patch for casinh.
author Markus Mützel <markus.muetzel@gmx.de>
date Mon, 08 Aug 2022 20:21:02 +0200
parents 6f3c099c0d38
children 1d99604afef8
files src/mingw-w64-complex-inverse-trig.patch
diffstat 1 files changed, 110 insertions(+), 38 deletions(-) [+]
line wrap: on
line diff
--- a/src/mingw-w64-complex-inverse-trig.patch	Sun Aug 07 20:09:23 2022 +0200
+++ b/src/mingw-w64-complex-inverse-trig.patch	Mon Aug 08 20:21:02 2022 +0200
@@ -120,28 +120,36 @@
 -- 
 2.35.3.windows.1
 
-From 98d5b204fe689ac35d348a1742723ec561a2c5f7 Mon Sep 17 00:00:00 2001
+From cddd2eb87270f368d18f79a1070a36e4dc5464eb 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.
+Date: Mon, 8 Aug 2022 18:01:32 +0200
+Subject: [PATCH 1/3] casinh: Use symmetries also for default case.
 
 ---
- mingw-w64-crt/complex/casinh.def.h | 46 ++++++++++++++++--------------
- 1 file changed, 25 insertions(+), 21 deletions(-)
+ mingw-w64-crt/complex/casinh.def.h | 47 +++++++++++++++++-------------
+ 1 file changed, 26 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
+index cce2d864e..40ff67579 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)
+@@ -47,6 +47,7 @@ __FLT_ABI(casinh) (__FLT_TYPE __complex__ z)
+ {
+   __complex__ __FLT_TYPE ret;
+   __complex__ __FLT_TYPE x;
++  __FLT_TYPE arz, aiz;
+   int r_class = fpclassify (__real__ z);
+   int i_class = fpclassify (__imag__ z);
+ 
+@@ -89,32 +90,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);
++  arz = __FLT_ABI(fabs) (__real__ z);
++  aiz = __FLT_ABI(fabs) (__imag__ z);
 +
 +  if (arz >= __FLT_CST(1.0)/__FLT_EPSILON
 +      || aiz >= __FLT_CST(1.0)/__FLT_EPSILON)
@@ -194,20 +202,10 @@
 -- 
 2.35.3.windows.1
 
-From 099f1cd2d4261573d584aa415f47742fa977c773 Mon Sep 17 00:00:00 2001
+From f309692d67674526ec9a87ecae527cb6c207cf93 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
+Date: Mon, 8 Aug 2022 18:02:40 +0200
+Subject: [PATCH 2/3] casinh: Use approximations for small real part and
  absolute imaginary part < 1
 
 ---
@@ -215,13 +213,27 @@
  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
+index 40ff67579..8d61b1399 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)
+@@ -104,6 +104,37 @@ __FLT_ABI(casinh) (__FLT_TYPE __complex__ z)
      ret = __FLT_ABI(clog) (x);
      __real__ ret += M_LN2;
    }
++  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 if (arz*arz <= __FLT_EPSILON && aiz < __FLT_CST(1.0))
 +  {
 +    /* Taylor series expansion around arz=0 for z + sqrt(z*z + 1):
@@ -239,23 +251,83 @@
 +    __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
 
+From b2ad7c96b759e8bc3af504018b70c91c8b0f135c Mon Sep 17 00:00:00 2001
+From: =?UTF-8?q?Markus=20M=C3=BCtzel?= <markus.muetzel@gmx.de>
+Date: Mon, 8 Aug 2022 18:03:49 +0200
+Subject: [PATCH 3/3] catanh: Use approximations for small real part
+
+---
+ mingw-w64-crt/complex/catanh.def.h | 37 +++++++++++++++++++++++++-----
+ 1 file changed, 31 insertions(+), 6 deletions(-)
+
+From 7589dfe1f5478a3ba7b93ea6dc28fc0ffd271f12 Mon Sep 17 00:00:00 2001
+From: =?UTF-8?q?Markus=20M=C3=BCtzel?= <markus.muetzel@gmx.de>
+Date: Mon, 8 Aug 2022 18:26:20 +0200
+Subject: [PATCH 3/3] catanh: Use approximations for small real part
+
+---
+ mingw-w64-crt/complex/catanh.def.h | 37 +++++++++++++++++++++++++-----
+ 1 file changed, 31 insertions(+), 6 deletions(-)
+
+diff --git a/mingw-w64-crt/complex/catanh.def.h b/mingw-w64-crt/complex/catanh.def.h
+index 68949d139..1d46d4dcd 100644
+--- a/mingw-w64-crt/complex/catanh.def.h
++++ b/mingw-w64-crt/complex/catanh.def.h
+@@ -75,17 +75,42 @@ __FLT_ABI(catanh) (__FLT_TYPE __complex__ z)
+   if (r_class == FP_ZERO && i_class == FP_ZERO)
+     return z;
+ 
++  /* catanh(z) = 1/2 * clog(1+z) - 1/2 * clog(1-z) = 1/2 * clog((1+z)/(1-z)) */
++
++  /* Use identity clog(c) = 1/2*log(|c|^2) + i*arg(c) to calculate real and
++  imaginary parts separately. */
++
++  /* real part */
++  /* |c|^2 = (Im(z)^2 + (1+Re(z))^2)/(Im(z)^2 + (1-Re(z))^2) */
+   i2 = __imag__ z * __imag__ z;
+ 
+-  n = __FLT_CST(1.0) + __real__ z;
+-  n = i2 + n * n;
++  if (__FLT_ABI(fabs) (__real__ z) <= __FLT_EPSILON)
++  {
++    /* |c|^2 = 1 + 4*Re(z)/(1+Im(z)^2) + O(Re(z)^2)  (Taylor series) */
++    __real__ ret = __FLT_CST(0.25) *
++      __FLT_ABI(log1p) (__FLT_CST(4.0)*(__real__ z) / (__FLT_CST(1.0) + i2));
++  }
++  else if ((__real__ z)*(__real__ z) <= __FLT_EPSILON)
++  {
++    /* |c|^2 = 1 + 4*Re(z)/(1+Im(z)^2) + 8*Re(z)^2/(1+Im(z)^2)^2 + O(Re(z)^3)  (Taylor series) */
++    d = __real__ z / (__FLT_CST(1.0) + i2);
++    __real__ ret = __FLT_CST(0.25) *
++      __FLT_ABI(log1p) (__FLT_CST(4.0) * d * (__FLT_CST(1.0) + __FLT_CST(2.0) * d));
++  }
++  else
++  {
++    n = __FLT_CST(1.0) + __real__ z;
++    n = i2 + n * n;
+ 
+-  d = __FLT_CST(1.0) - __real__ z;
+-  d = i2 + d * d;
++    d = __FLT_CST(1.0) - __real__ z;
++    d = i2 + d * d;
+ 
+-  __real__ ret = __FLT_CST(0.25) * (__FLT_ABI(log) (n) - __FLT_ABI(log) (d));
++    __real__ ret = __FLT_CST(0.25) * (__FLT_ABI(log) (n) - __FLT_ABI(log) (d));
++  }
+ 
+-  d = 1 - __real__ z * __real__ z - i2;
++  /* imaginary part */
++  /* z = (1 - Re(z)^2 - Im(z)^2 + 2i * Im(z) / ((1-Re(z))^2 + Im(z)^2) */
++  d = __FLT_CST(1.0) - __real__ z * __real__ z - i2;
+ 
+   __imag__ ret = __FLT_CST(0.5) * __FLT_ABI(atan2) (__FLT_CST(2.0) * __imag__ z, d);
+ 
+-- 
+2.35.3.windows.1
+