# HG changeset patch # User Markus Mützel # Date 1659895763 -7200 # Node ID 6f3c099c0d38e191e232f75ff6ef6e169e9262f1 # Parent 947f8511e9d1dd51ffead86cd7b7e1df8e93b710 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. diff -r 947f8511e9d1 -r 6f3c099c0d38 src/mingw-w64-complex-inverse-trig.patch --- a/src/mingw-w64-complex-inverse-trig.patch Fri Aug 05 19:40:46 2022 -0400 +++ b/src/mingw-w64-complex-inverse-trig.patch Sun Aug 07 20:09:23 2022 +0200 @@ -120,3 +120,142 @@ -- 2.35.3.windows.1 +From 98d5b204fe689ac35d348a1742723ec561a2c5f7 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. + +--- + 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?= +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 + 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 +