Mercurial > mxe-octave
view src/mingw-w64-complex-inverse-trig.patch @ 6503:bf4a0944be3a
* src/mpg123.mk: update to v1.31.1
author | John Donoghue <john.donoghue@ieee.org> |
---|---|
date | Thu, 03 Nov 2022 16:24:05 -0400 |
parents | dbc29160e2b2 |
children |
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 (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 (signbit (__real__ ret)) ret = -ret; return ret; -- 2.35.3.windows.1 From cddd2eb87270f368d18f79a1070a36e4dc5464eb 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:01:32 +0200 Subject: [PATCH 1/3] casinh: Use symmetries also for default case. --- 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..40ff67579 100644 --- a/mingw-w64-crt/complex/casinh.def.h +++ b/mingw-w64-crt/complex/casinh.def.h @@ -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. */ + 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) { /* 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 f309692d67674526ec9a87ecae527cb6c207cf93 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:02:40 +0200 Subject: [PATCH 2/3] 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 40ff67579..8d61b1399 100644 --- a/mingw-w64-crt/complex/casinh.def.h +++ b/mingw-w64-crt/complex/casinh.def.h @@ -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): + 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 { __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