Mercurial > mxe-octave
comparison src/mingw-w64-complex-inverse-trig.patch @ 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 | dbc29160e2b2 |
comparison
equal
deleted
inserted
replaced
6356:6f3c099c0d38 | 6357:dc5ad8056086 |
---|---|
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 | 123 From cddd2eb87270f368d18f79a1070a36e4dc5464eb Mon Sep 17 00:00:00 2001 |
124 From: =?UTF-8?q?Markus=20M=C3=BCtzel?= <markus.muetzel@gmx.de> | 124 From: =?UTF-8?q?Markus=20M=C3=BCtzel?= <markus.muetzel@gmx.de> |
125 Date: Sun, 7 Aug 2022 13:21:29 +0200 | 125 Date: Mon, 8 Aug 2022 18:01:32 +0200 |
126 Subject: [PATCH 1/2] casinh: Use symmetries also for default case. | 126 Subject: [PATCH 1/3] casinh: Use symmetries also for default case. |
127 | 127 |
128 --- | 128 --- |
129 mingw-w64-crt/complex/casinh.def.h | 46 ++++++++++++++++-------------- | 129 mingw-w64-crt/complex/casinh.def.h | 47 +++++++++++++++++------------- |
130 1 file changed, 25 insertions(+), 21 deletions(-) | 130 1 file changed, 26 insertions(+), 21 deletions(-) |
131 | 131 |
132 diff --git a/mingw-w64-crt/complex/casinh.def.h b/mingw-w64-crt/complex/casinh.def.h | 132 diff --git a/mingw-w64-crt/complex/casinh.def.h b/mingw-w64-crt/complex/casinh.def.h |
133 index cce2d864e..d309e6392 100644 | 133 index cce2d864e..40ff67579 100644 |
134 --- a/mingw-w64-crt/complex/casinh.def.h | 134 --- a/mingw-w64-crt/complex/casinh.def.h |
135 +++ b/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) | 136 @@ -47,6 +47,7 @@ __FLT_ABI(casinh) (__FLT_TYPE __complex__ z) |
137 { | |
138 __complex__ __FLT_TYPE ret; | |
139 __complex__ __FLT_TYPE x; | |
140 + __FLT_TYPE arz, aiz; | |
141 int r_class = fpclassify (__real__ z); | |
142 int i_class = fpclassify (__imag__ z); | |
143 | |
144 @@ -89,32 +90,36 @@ __FLT_ABI(casinh) (__FLT_TYPE __complex__ z) | |
137 | 145 |
138 /* casinh(z) = log(z + sqrt(z*z + 1)) */ | 146 /* casinh(z) = log(z + sqrt(z*z + 1)) */ |
139 | 147 |
140 - if (__FLT_ABI(fabs) (__real__ z) >= __FLT_CST(1.0)/__FLT_EPSILON | 148 - if (__FLT_ABI(fabs) (__real__ z) >= __FLT_CST(1.0)/__FLT_EPSILON |
141 - || __FLT_ABI(fabs) (__imag__ z) >= __FLT_CST(1.0)/__FLT_EPSILON) | 149 - || __FLT_ABI(fabs) (__imag__ z) >= __FLT_CST(1.0)/__FLT_EPSILON) |
142 + /* Use symmetries to perform the calculation in the first quadrant. */ | 150 + /* Use symmetries to perform the calculation in the first quadrant. */ |
143 + __FLT_TYPE arz = __FLT_ABI(fabs) (__real__ z); | 151 + arz = __FLT_ABI(fabs) (__real__ z); |
144 + __FLT_TYPE aiz = __FLT_ABI(fabs) (__imag__ z); | 152 + aiz = __FLT_ABI(fabs) (__imag__ z); |
145 + | 153 + |
146 + if (arz >= __FLT_CST(1.0)/__FLT_EPSILON | 154 + if (arz >= __FLT_CST(1.0)/__FLT_EPSILON |
147 + || aiz >= __FLT_CST(1.0)/__FLT_EPSILON) | 155 + || aiz >= __FLT_CST(1.0)/__FLT_EPSILON) |
148 { | 156 { |
149 /* For large z, z + sqrt(z*z + 1) is approximately 2*z. | 157 /* For large z, z + sqrt(z*z + 1) is approximately 2*z. |
192 + return ret; | 200 + return ret; |
193 } | 201 } |
194 -- | 202 -- |
195 2.35.3.windows.1 | 203 2.35.3.windows.1 |
196 | 204 |
197 From 099f1cd2d4261573d584aa415f47742fa977c773 Mon Sep 17 00:00:00 2001 | 205 From f309692d67674526ec9a87ecae527cb6c207cf93 Mon Sep 17 00:00:00 2001 |
198 From: =?UTF-8?q?Markus=20M=C3=BCtzel?= <markus.muetzel@gmx.de> | 206 From: =?UTF-8?q?Markus=20M=C3=BCtzel?= <markus.muetzel@gmx.de> |
199 Date: Sun, 7 Aug 2022 13:23:02 +0200 | 207 Date: Mon, 8 Aug 2022 18:02:40 +0200 |
200 Subject: [PATCH 2/2] casinh: Use approximation for small real part and | 208 Subject: [PATCH 2/3] casinh: Use approximations for small real part and |
201 absolute imaginary part < 1 | 209 absolute imaginary part < 1 |
202 | 210 |
203 --- | 211 --- |
204 mingw-w64-crt/complex/casinh.def.h | 31 ++++++++++++++++++++++++++++++ | 212 mingw-w64-crt/complex/casinh.def.h | 31 ++++++++++++++++++++++++++++++ |
205 1 file changed, 31 insertions(+) | 213 1 file changed, 31 insertions(+) |
206 | 214 |
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 | 215 diff --git a/mingw-w64-crt/complex/casinh.def.h b/mingw-w64-crt/complex/casinh.def.h |
218 index d309e6392..a5447bfcd 100644 | 216 index 40ff67579..8d61b1399 100644 |
219 --- a/mingw-w64-crt/complex/casinh.def.h | 217 --- a/mingw-w64-crt/complex/casinh.def.h |
220 +++ b/mingw-w64-crt/complex/casinh.def.h | 218 +++ b/mingw-w64-crt/complex/casinh.def.h |
221 @@ -103,6 +103,37 @@ __FLT_ABI(casinh) (__FLT_TYPE __complex__ z) | 219 @@ -104,6 +104,37 @@ __FLT_ABI(casinh) (__FLT_TYPE __complex__ z) |
222 ret = __FLT_ABI(clog) (x); | 220 ret = __FLT_ABI(clog) (x); |
223 __real__ ret += M_LN2; | 221 __real__ ret += M_LN2; |
224 } | 222 } |
223 + else if (arz <= __FLT_EPSILON && aiz < __FLT_CST(1.0)) | |
224 + { | |
225 + /* Taylor series expansion around arz=0 for z + sqrt(z*z + 1): | |
226 + c = arz + sqrt(1-aiz^2) + i*(aiz + arz*aiz / sqrt(1-aiz^2)) + O(arz^2) | |
227 + Identity: clog(c) = log(|c|) + i*arg(c) | |
228 + For real part of result: | |
229 + |c| = 1 + arz / sqrt(1-aiz^2) + O(arz^2) (Taylor series expansion) | |
230 + For imaginary part of result: | |
231 + c = (arz + sqrt(1-aiz^2))/sqrt(1-aiz^2) * (sqrt(1-aiz^2) + i*aiz) + O(arz^6) | |
232 + */ | |
233 + __FLT_TYPE s1maiz2 = __FLT_ABI(sqrt) ((__FLT_CST(1.0)+aiz)*(__FLT_CST(1.0)-aiz)); | |
234 + __real__ ret = __FLT_ABI(log1p) (arz / s1maiz2); | |
235 + __imag__ ret = __FLT_ABI(atan2) (aiz, s1maiz2); | |
236 + } | |
225 + else if (arz*arz <= __FLT_EPSILON && aiz < __FLT_CST(1.0)) | 237 + else if (arz*arz <= __FLT_EPSILON && aiz < __FLT_CST(1.0)) |
226 + { | 238 + { |
227 + /* Taylor series expansion around arz=0 for z + sqrt(z*z + 1): | 239 + /* 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) | 240 + 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) | 241 + Identity: clog(c) = log(|c|) + i*arg(c) |
237 + __FLT_TYPE arz2red = arz * arz / __FLT_CST(2.0) / s1maiz2; | 249 + __FLT_TYPE arz2red = arz * arz / __FLT_CST(2.0) / s1maiz2; |
238 + __real__ ret = __FLT_ABI(log1p) ((arz + arz2red) / s1maiz2); | 250 + __real__ ret = __FLT_ABI(log1p) ((arz + arz2red) / s1maiz2); |
239 + __imag__ ret = __FLT_ABI(atan2) (aiz * (s1maiz2 + arz), | 251 + __imag__ ret = __FLT_ABI(atan2) (aiz * (s1maiz2 + arz), |
240 + onemaiz2 + arz*s1maiz2 + arz2red); | 252 + onemaiz2 + arz*s1maiz2 + arz2red); |
241 + } | 253 + } |
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 | 254 else |
257 { | 255 { |
258 __real__ x = (arz - aiz) * (arz + aiz) + __FLT_CST(1.0); | 256 __real__ x = (arz - aiz) * (arz + aiz) + __FLT_CST(1.0); |
259 -- | 257 -- |
260 2.35.3.windows.1 | 258 2.35.3.windows.1 |
261 | 259 |
260 From b2ad7c96b759e8bc3af504018b70c91c8b0f135c Mon Sep 17 00:00:00 2001 | |
261 From: =?UTF-8?q?Markus=20M=C3=BCtzel?= <markus.muetzel@gmx.de> | |
262 Date: Mon, 8 Aug 2022 18:03:49 +0200 | |
263 Subject: [PATCH 3/3] catanh: Use approximations for small real part | |
264 | |
265 --- | |
266 mingw-w64-crt/complex/catanh.def.h | 37 +++++++++++++++++++++++++----- | |
267 1 file changed, 31 insertions(+), 6 deletions(-) | |
268 | |
269 From 7589dfe1f5478a3ba7b93ea6dc28fc0ffd271f12 Mon Sep 17 00:00:00 2001 | |
270 From: =?UTF-8?q?Markus=20M=C3=BCtzel?= <markus.muetzel@gmx.de> | |
271 Date: Mon, 8 Aug 2022 18:26:20 +0200 | |
272 Subject: [PATCH 3/3] catanh: Use approximations for small real part | |
273 | |
274 --- | |
275 mingw-w64-crt/complex/catanh.def.h | 37 +++++++++++++++++++++++++----- | |
276 1 file changed, 31 insertions(+), 6 deletions(-) | |
277 | |
278 diff --git a/mingw-w64-crt/complex/catanh.def.h b/mingw-w64-crt/complex/catanh.def.h | |
279 index 68949d139..1d46d4dcd 100644 | |
280 --- a/mingw-w64-crt/complex/catanh.def.h | |
281 +++ b/mingw-w64-crt/complex/catanh.def.h | |
282 @@ -75,17 +75,42 @@ __FLT_ABI(catanh) (__FLT_TYPE __complex__ z) | |
283 if (r_class == FP_ZERO && i_class == FP_ZERO) | |
284 return z; | |
285 | |
286 + /* catanh(z) = 1/2 * clog(1+z) - 1/2 * clog(1-z) = 1/2 * clog((1+z)/(1-z)) */ | |
287 + | |
288 + /* Use identity clog(c) = 1/2*log(|c|^2) + i*arg(c) to calculate real and | |
289 + imaginary parts separately. */ | |
290 + | |
291 + /* real part */ | |
292 + /* |c|^2 = (Im(z)^2 + (1+Re(z))^2)/(Im(z)^2 + (1-Re(z))^2) */ | |
293 i2 = __imag__ z * __imag__ z; | |
294 | |
295 - n = __FLT_CST(1.0) + __real__ z; | |
296 - n = i2 + n * n; | |
297 + if (__FLT_ABI(fabs) (__real__ z) <= __FLT_EPSILON) | |
298 + { | |
299 + /* |c|^2 = 1 + 4*Re(z)/(1+Im(z)^2) + O(Re(z)^2) (Taylor series) */ | |
300 + __real__ ret = __FLT_CST(0.25) * | |
301 + __FLT_ABI(log1p) (__FLT_CST(4.0)*(__real__ z) / (__FLT_CST(1.0) + i2)); | |
302 + } | |
303 + else if ((__real__ z)*(__real__ z) <= __FLT_EPSILON) | |
304 + { | |
305 + /* |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) */ | |
306 + d = __real__ z / (__FLT_CST(1.0) + i2); | |
307 + __real__ ret = __FLT_CST(0.25) * | |
308 + __FLT_ABI(log1p) (__FLT_CST(4.0) * d * (__FLT_CST(1.0) + __FLT_CST(2.0) * d)); | |
309 + } | |
310 + else | |
311 + { | |
312 + n = __FLT_CST(1.0) + __real__ z; | |
313 + n = i2 + n * n; | |
314 | |
315 - d = __FLT_CST(1.0) - __real__ z; | |
316 - d = i2 + d * d; | |
317 + d = __FLT_CST(1.0) - __real__ z; | |
318 + d = i2 + d * d; | |
319 | |
320 - __real__ ret = __FLT_CST(0.25) * (__FLT_ABI(log) (n) - __FLT_ABI(log) (d)); | |
321 + __real__ ret = __FLT_CST(0.25) * (__FLT_ABI(log) (n) - __FLT_ABI(log) (d)); | |
322 + } | |
323 | |
324 - d = 1 - __real__ z * __real__ z - i2; | |
325 + /* imaginary part */ | |
326 + /* z = (1 - Re(z)^2 - Im(z)^2 + 2i * Im(z) / ((1-Re(z))^2 + Im(z)^2) */ | |
327 + d = __FLT_CST(1.0) - __real__ z * __real__ z - i2; | |
328 | |
329 __imag__ ret = __FLT_CST(0.5) * __FLT_ABI(atan2) (__FLT_CST(2.0) * __imag__ z, d); | |
330 | |
331 -- | |
332 2.35.3.windows.1 | |
333 |