Mercurial > mxe-octave
comparison src/mingw-w64-complex-inverse-trig.patch @ 6356:6f3c099c0d38
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 |
comparison
equal
deleted
inserted
replaced
6355:947f8511e9d1 | 6356:6f3c099c0d38 |
---|---|
118 | 118 |
119 return ret; | 119 return ret; |
120 -- | 120 -- |
121 2.35.3.windows.1 | 121 2.35.3.windows.1 |
122 | 122 |
123 From 98d5b204fe689ac35d348a1742723ec561a2c5f7 Mon Sep 17 00:00:00 2001 | |
124 From: =?UTF-8?q?Markus=20M=C3=BCtzel?= <markus.muetzel@gmx.de> | |
125 Date: Sun, 7 Aug 2022 13:21:29 +0200 | |
126 Subject: [PATCH 1/2] casinh: Use symmetries also for default case. | |
127 | |
128 --- | |
129 mingw-w64-crt/complex/casinh.def.h | 46 ++++++++++++++++-------------- | |
130 1 file changed, 25 insertions(+), 21 deletions(-) | |
131 | |
132 diff --git a/mingw-w64-crt/complex/casinh.def.h b/mingw-w64-crt/complex/casinh.def.h | |
133 index cce2d864e..d309e6392 100644 | |
134 --- a/mingw-w64-crt/complex/casinh.def.h | |
135 +++ b/mingw-w64-crt/complex/casinh.def.h | |
136 @@ -89,32 +89,36 @@ __FLT_ABI(casinh) (__FLT_TYPE __complex__ z) | |
137 | |
138 /* casinh(z) = log(z + sqrt(z*z + 1)) */ | |
139 | |
140 - if (__FLT_ABI(fabs) (__real__ z) >= __FLT_CST(1.0)/__FLT_EPSILON | |
141 - || __FLT_ABI(fabs) (__imag__ z) >= __FLT_CST(1.0)/__FLT_EPSILON) | |
142 + /* Use symmetries to perform the calculation in the first quadrant. */ | |
143 + __FLT_TYPE arz = __FLT_ABI(fabs) (__real__ z); | |
144 + __FLT_TYPE aiz = __FLT_ABI(fabs) (__imag__ z); | |
145 + | |
146 + if (arz >= __FLT_CST(1.0)/__FLT_EPSILON | |
147 + || aiz >= __FLT_CST(1.0)/__FLT_EPSILON) | |
148 { | |
149 /* For large z, z + sqrt(z*z + 1) is approximately 2*z. | |
150 - Use that approximation to avoid overflow when squaring. | |
151 - Additionally, use symmetries to perform the calculation in the first | |
152 - quadrant. */ | |
153 - __real__ x = __FLT_ABI(fabs) (__real__ z); | |
154 - __imag__ x = __FLT_ABI(fabs) (__imag__ z); | |
155 - x = __FLT_ABI(clog) (x); | |
156 - __real__ x += M_LN2; | |
157 - | |
158 - /* adjust signs for input quadrant */ | |
159 - __real__ ret = __FLT_ABI(copysign) (__real__ x, __real__ z); | |
160 - __imag__ ret = __FLT_ABI(copysign) (__imag__ x, __imag__ z); | |
161 - | |
162 - return ret; | |
163 + Use that approximation to avoid overflow when squaring. */ | |
164 + __real__ x = arz; | |
165 + __imag__ x = aiz; | |
166 + ret = __FLT_ABI(clog) (x); | |
167 + __real__ ret += M_LN2; | |
168 } | |
169 + else | |
170 + { | |
171 + __real__ x = (arz - aiz) * (arz + aiz) + __FLT_CST(1.0); | |
172 + __imag__ x = __FLT_CST(2.0) * arz * aiz; | |
173 + | |
174 + x = __FLT_ABI(csqrt) (x); | |
175 | |
176 - __real__ x = (__real__ z - __imag__ z) * (__real__ z + __imag__ z) + __FLT_CST(1.0); | |
177 - __imag__ x = __FLT_CST(2.0) * __real__ z * __imag__ z; | |
178 + __real__ x += arz; | |
179 + __imag__ x += aiz; | |
180 | |
181 - x = __FLT_ABI(csqrt) (x); | |
182 + ret = __FLT_ABI(clog) (x); | |
183 + } | |
184 | |
185 - __real__ x += __real__ z; | |
186 - __imag__ x += __imag__ z; | |
187 + /* adjust signs for input quadrant */ | |
188 + __real__ ret = __FLT_ABI(copysign) (__real__ ret, __real__ z); | |
189 + __imag__ ret = __FLT_ABI(copysign) (__imag__ ret, __imag__ z); | |
190 | |
191 - return __FLT_ABI(clog) (x); | |
192 + return ret; | |
193 } | |
194 -- | |
195 2.35.3.windows.1 | |
196 | |
197 From 099f1cd2d4261573d584aa415f47742fa977c773 Mon Sep 17 00:00:00 2001 | |
198 From: =?UTF-8?q?Markus=20M=C3=BCtzel?= <markus.muetzel@gmx.de> | |
199 Date: Sun, 7 Aug 2022 13:23:02 +0200 | |
200 Subject: [PATCH 2/2] casinh: Use approximation for small real part and | |
201 absolute imaginary part < 1 | |
202 | |
203 --- | |
204 mingw-w64-crt/complex/casinh.def.h | 31 ++++++++++++++++++++++++++++++ | |
205 1 file changed, 31 insertions(+) | |
206 | |
207 From 575a070e3fb360fced10a458e72a298d3071c42e Mon Sep 17 00:00:00 2001 | |
208 From: =?UTF-8?q?Markus=20M=C3=BCtzel?= <markus.muetzel@gmx.de> | |
209 Date: Sun, 7 Aug 2022 19:05:21 +0200 | |
210 Subject: [PATCH 2/2] casinh: Use approximations for small real part and | |
211 absolute imaginary part < 1 | |
212 | |
213 --- | |
214 mingw-w64-crt/complex/casinh.def.h | 31 ++++++++++++++++++++++++++++++ | |
215 1 file changed, 31 insertions(+) | |
216 | |
217 diff --git a/mingw-w64-crt/complex/casinh.def.h b/mingw-w64-crt/complex/casinh.def.h | |
218 index d309e6392..a5447bfcd 100644 | |
219 --- a/mingw-w64-crt/complex/casinh.def.h | |
220 +++ b/mingw-w64-crt/complex/casinh.def.h | |
221 @@ -103,6 +103,37 @@ __FLT_ABI(casinh) (__FLT_TYPE __complex__ z) | |
222 ret = __FLT_ABI(clog) (x); | |
223 __real__ ret += M_LN2; | |
224 } | |
225 + else if (arz*arz <= __FLT_EPSILON && aiz < __FLT_CST(1.0)) | |
226 + { | |
227 + /* Taylor series expansion around arz=0 for z + sqrt(z*z + 1): | |
228 + 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) | |
229 + Identity: clog(c) = log(|c|) + i*arg(c) | |
230 + For real part of result: | |
231 + |c| = 1 + arz / sqrt(1-aiz^2) + arz^2/(2*(1-aiz^2)) + O(arz^3) (Taylor series expansion) | |
232 + For imaginary part of result: | |
233 + 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) | |
234 + */ | |
235 + __FLT_TYPE onemaiz2 = (__FLT_CST(1.0)+aiz)*(__FLT_CST(1.0)-aiz); | |
236 + __FLT_TYPE s1maiz2 = __FLT_ABI(sqrt) (onemaiz2); | |
237 + __FLT_TYPE arz2red = arz * arz / __FLT_CST(2.0) / s1maiz2; | |
238 + __real__ ret = __FLT_ABI(log1p) ((arz + arz2red) / s1maiz2); | |
239 + __imag__ ret = __FLT_ABI(atan2) (aiz * (s1maiz2 + arz), | |
240 + onemaiz2 + arz*s1maiz2 + arz2red); | |
241 + } | |
242 + else if (arz <= __FLT_EPSILON && aiz < __FLT_CST(1.0)) | |
243 + { | |
244 + /* Taylor series expansion around arz=0 for z + sqrt(z*z + 1): | |
245 + c = arz + sqrt(1-aiz^2) + i*(aiz + arz*aiz / sqrt(1-aiz^2)) + O(arz^2) | |
246 + Identity: clog(c) = log(|c|) + i*arg(c) | |
247 + For real part of result: | |
248 + |c| = 1 + arz / sqrt(1-aiz^2) + O(arz^2) (Taylor series expansion) | |
249 + For imaginary part of result: | |
250 + c = (arz + sqrt(1-aiz^2))/sqrt(1-aiz^2) * (sqrt(1-aiz^2) + i*aiz) + O(arz^6) | |
251 + */ | |
252 + __FLT_TYPE s1maiz2 = __FLT_ABI(sqrt) ((__FLT_CST(1.0)+aiz)*(__FLT_CST(1.0)-aiz)); | |
253 + __real__ ret = __FLT_ABI(log1p) (arz / s1maiz2); | |
254 + __imag__ ret = __FLT_ABI(atan2) (aiz, s1maiz2); | |
255 + } | |
256 else | |
257 { | |
258 __real__ x = (arz - aiz) * (arz + aiz) + __FLT_CST(1.0); | |
259 -- | |
260 2.35.3.windows.1 | |
261 |