# HG changeset patch # User Markus Mützel # Date 1659982862 -7200 # Node ID dc5ad80560861dac66d8cc1d3389025f1604c39b # Parent 6f3c099c0d38e191e232f75ff6ef6e169e9262f1 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. diff -r 6f3c099c0d38 -r dc5ad8056086 src/mingw-w64-complex-inverse-trig.patch --- 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?= -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?= -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?= -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?= +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?= +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 +