Oops. Resending patches in .txt files.
I’ve also reattached the test file, to which I added a small amount of
additional testing.
Thanks, David.
From: JonY via Mingw-w64-public<mailto:[email protected]>
Sent: 03 October 2021 04:12
To:
[email protected]<mailto:[email protected]>
Cc: JonY<mailto:[email protected]>
Subject: Re: [Mingw-w64-public] Patches to fix #515 and 916.
On 10/2/21 5:27 PM, David James wrote:
>
> Hello – I’ve coded some fixed for bugs #515 Incorrect sign output from
> asinh<https://sourceforge.net/p/mingw-w64/bugs/515/> (which also impacted
> atanh) and #916 asinh incorrect for large
> values<https://sourceforge.net/p/mingw-w64/bugs/916/> and attached some patch
> files. I’m hoping to get these fixes into GHC Haskell (I believe they need to
> get to msys2 first).
>
> Please could you review / comment?
>
> I’ve also included a file I used for testing the fixes. (I don’t know if this
> could/should be integrated into any CI testing?)
>
> Please let me know if there are any questions, etc. (It’s the first time I’ve
> created patch files, so please let me know if I’ve done something wrong).
>
> Thanks!
> David.
>
Looks like SF ate your patches, please resend with the .txt extension,
it may help.
_______________________________________________
Mingw-w64-public mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/mingw-w64-public
From f03ed1372d031e529ae9a2cdaad4db7ea66c95af Mon Sep 17 00:00:00 2001
From: davjam <[email protected]>
Date: Sat, 2 Oct 2021 17:33:26 +0100
Subject: [PATCH 1/3] Fixes asinh,f,l for bugs #515 and #916
ie. fix sign of asinh(-0) and ensure results returned for large values.
---
mingw-w64-crt/math/x86/asinh.c | 55 ++++++++++++++++++++++++++++-----
mingw-w64-crt/math/x86/asinhf.c | 17 +++++-----
mingw-w64-crt/math/x86/asinhl.c | 16 +++++-----
3 files changed, 65 insertions(+), 23 deletions(-)
diff --git a/mingw-w64-crt/math/x86/asinh.c b/mingw-w64-crt/math/x86/asinh.c
index 826cc6dee..002aedee8 100644
--- a/mingw-w64-crt/math/x86/asinh.c
+++ b/mingw-w64-crt/math/x86/asinh.c
@@ -5,6 +5,7 @@
*/
#include <math.h>
#include <errno.h>
+#include <float.h>
#include "fastmath.h"
/* asinh(x) = copysign(log(fabs(x) + sqrt(x * x + 1.0)), x) */
@@ -21,13 +22,53 @@ double asinh(double x)
return x;
#endif
- /* Use log1p to avoid cancellation with small x. Put
- x * x in denom, so overflow is harmless.
- asinh(x) = log1p (x + sqrt (x * x + 1.0) - 1.0)
- = log1p (x + x * x / (sqrt (x * x + 1.0) + 1.0)) */
+ /* NB the previous formula
+ z = __fast_log1p (z + z * z / (__fast_sqrt (z * z + 1.0) + 1.0));
+ was defective in two ways:
+ 1: It ommitted required brackets:
+ z = __fast_log1p (z + z * (z / (__fast_sqrt (z * z + 1.0) + 1.0)));
+ ^ ^
+ so would still overflow for large z.
+ 2: Even with the brackets, it still degraded quickly for large z
+ (where z*z+1 == z*z).
+ e.g. asinh (sinh 356.0)) gave 355.30685281944005
+ */
- z = __fast_log1p (z + z * z / (__fast_sqrt (z * z + 1.0) + 1.0));
+ const double asinhCutover = pow(2,DBL_MAX_EXP/2); // 1.3407807929943e+154
- return ( x > 0.0 ? z : -z);
+ if (z < asinhCutover)
+ /* After excluding large values, the rearranged formula gives better results
+ the original formula log(z + sqrt(z * z + 1.0)) for very small z.
+ e.g. rearranged asinh(sinh 2e-301)) = 2e-301
+ original asinh(sinh 2e-301)) = 0.
+ asinh(z) = log (z + sqrt (z * z + 1.0))
+ = log1p (z + sqrt (z * z + 1.0) - 1.0)
+ = log1p (z + (sqrt (z * z + 1.0) - 1.0)
+ * (sqrt (z * z + 1.0) + 1.0)
+ / (sqrt (z * z + 1.0) + 1.0))
+ = log1p (z + ((z * z + 1.0) - 1.0)
+ / (sqrt (z * z + 1.0) + 1.0))
+ = log1p (z + z * z / (sqrt (z * z + 1.0) + 1.0))
+ */
+ z = __fast_log1p (z + z * (z / (__fast_sqrt (z * z + 1.0) + 1.0)));
+ else
+ /* above this, z*z+1 == z*z, so we can simplify
+ (and avoid z*z being infinity).
+ asinh(z) = log (z + sqrt (z * z + 1.0))
+ = log (z + sqrt (z * z ))
+ = log (2 * z)
+ = log 2 + log z
+ Choosing asinhCutover is a little tricky.
+ We'd like something that's based on the nature of
+ the numeric type (DBL_MAX_EXP, etc).
+ If c = asinhCutover, then we need:
+ (1) c*c == c*c + 1
+ (2) log (2*c) = log 2 + log c.
+ For float:
+ 9.490626562425156e7 is the smallest value that
+ achieves (1), but it fails (2). (It only just fails,
+ but enough to make the function erroneously non-monotonic).
+ */
+ z = __fast_log(2) + __fast_log(z);
+ return copysign(z, x); //ensure 0.0 -> 0.0 and -0.0 -> -0.0.
}
-
diff --git a/mingw-w64-crt/math/x86/asinhf.c b/mingw-w64-crt/math/x86/asinhf.c
index fae785abe..f27f84314 100644
--- a/mingw-w64-crt/math/x86/asinhf.c
+++ b/mingw-w64-crt/math/x86/asinhf.c
@@ -5,6 +5,7 @@
*/
#include <math.h>
#include <errno.h>
+#include <float.h>
#include "fastmath.h"
/* asinh(x) = copysign(log(fabs(x) + sqrt(x * x + 1.0)), x) */
@@ -21,13 +22,13 @@ float asinhf(float x)
return x;
#endif
+ /* See commentary in asinh */
+ const float asinhCutover = pow(2,FLT_MAX_EXP/2);
- /* Use log1p to avoid cancellation with small x. Put
- x * x in denom, so overflow is harmless.
- asinh(x) = log1p (x + sqrt (x * x + 1.0) - 1.0)
- = log1p (x + x * x / (sqrt (x * x + 1.0) + 1.0)) */
-
- z = __fast_log1p (z + z * z / (__fast_sqrt (z * z + 1.0) + 1.0));
-
- return ( x > 0.0 ? z : -z);
+ if (z < asinhCutover)
+ z = __fast_log1p (z + z * (z / (__fast_sqrt (z * z + 1.0) + 1.0)));
+ //z = __fast_log(z + __fast_sqrt(z * z + 1.0));
+ else
+ z = __fast_log(2) + __fast_log(z);
+ return copysignf(z, x); //ensure 0.0 -> 0.0 and -0.0 -> -0.0.
}
diff --git a/mingw-w64-crt/math/x86/asinhl.c b/mingw-w64-crt/math/x86/asinhl.c
index bb2ca97b2..7c0a64e2a 100644
--- a/mingw-w64-crt/math/x86/asinhl.c
+++ b/mingw-w64-crt/math/x86/asinhl.c
@@ -5,6 +5,7 @@
*/
#include <math.h>
#include <errno.h>
+#include <float.h>
#include "fastmath.h"
/* asinh(x) = copysign(log(fabs(x) + sqrt(x * x + 1.0)), x) */
@@ -13,7 +14,6 @@ long double asinhl(long double x)
long double z;
if (!isfinite (x))
return x;
-
z = fabsl (x);
/* Avoid setting FPU underflow exception flag in x * x. */
@@ -22,12 +22,12 @@ long double asinhl(long double x)
return x;
#endif
- /* Use log1p to avoid cancellation with small x. Put
- x * x in denom, so overflow is harmless.
- asinh(x) = log1p (x + sqrt (x * x + 1.0) - 1.0)
- = log1p (x + x * x / (sqrt (x * x + 1.0) + 1.0)) */
-
- z = __fast_log1pl (z + z * z / (__fast_sqrtl (z * z + 1.0L) + 1.0L));
+ /* See commentary in asinh */
+ const long double asinhCutover = powl(2,LDBL_MAX_EXP/2);
- return ( x > 0.0 ? z : -z);
+ if (z < asinhCutover)
+ z = __fast_log1pl (z + z * (z / (__fast_sqrtl (z * z + 1.0) + 1.0)));
+ else
+ z = __fast_logl(2) + __fast_logl(z);
+ return copysignl(z, x); //ensure 0.0 -> 0.0 and -0.0 -> -0.0.
}
--
2.25.1.windows.1
From 1ac64e22b5e30ab543b37f142cef517a624c576d Mon Sep 17 00:00:00 2001
From: davjam <[email protected]>
Date: Sat, 2 Oct 2021 17:34:13 +0100
Subject: [PATCH 2/3] fix atanh,f,l for bug #515
atanh(-0) now gives -0.
---
mingw-w64-crt/math/x86/atanh.c | 2 +-
mingw-w64-crt/math/x86/atanhf.c | 2 +-
mingw-w64-crt/math/x86/atanhl.c | 2 +-
3 files changed, 3 insertions(+), 3 deletions(-)
diff --git a/mingw-w64-crt/math/x86/atanh.c b/mingw-w64-crt/math/x86/atanh.c
index 429bfeae7..b3aa76fe9 100644
--- a/mingw-w64-crt/math/x86/atanh.c
+++ b/mingw-w64-crt/math/x86/atanh.c
@@ -32,5 +32,5 @@ double atanh(double x)
= 0.5 * log1p ((1.0 + x - 1.0 + x) /(1.0 - x))
= 0.5 * log1p ((2.0 * x ) / (1.0 - x)) */
z = 0.5 * __fast_log1p ((z + z) / (1.0 - z));
- return x >= 0 ? z : -z;
+ return copysign(z, x); //ensure 0.0 -> 0.0 and -0.0 -> -0.0.
}
diff --git a/mingw-w64-crt/math/x86/atanhf.c b/mingw-w64-crt/math/x86/atanhf.c
index 96d1e9b22..9adc6902d 100644
--- a/mingw-w64-crt/math/x86/atanhf.c
+++ b/mingw-w64-crt/math/x86/atanhf.c
@@ -31,5 +31,5 @@ float atanhf (float x)
= 0.5 * log1p ((1.0 + x - 1.0 + x) /(1.0 - x))
= 0.5 * log1p ((2.0 * x ) / (1.0 - x)) */
z = 0.5 * __fast_log1p ((z + z) / (1.0 - z));
- return x >= 0 ? z : -z;
+ return copysignf(z, x); //ensure 0.0 -> 0.0 and -0.0 -> -0.0.
}
diff --git a/mingw-w64-crt/math/x86/atanhl.c b/mingw-w64-crt/math/x86/atanhl.c
index 59eb1bd91..76897466e 100644
--- a/mingw-w64-crt/math/x86/atanhl.c
+++ b/mingw-w64-crt/math/x86/atanhl.c
@@ -30,5 +30,5 @@ long double atanhl (long double x)
= 0.5 * log1p ((1.0 + x - 1.0 + x) /(1.0 - x))
= 0.5 * log1p ((2.0 * x ) / (1.0 - x)) */
z = 0.5L * __fast_log1pl ((z + z) / (1.0L - z));
- return x >= 0 ? z : -z;
+ return copysignl(z, x); //ensure 0.0 -> 0.0 and -0.0 -> -0.0.
}
--
2.25.1.windows.1
#include <stdio.h>
#include <string.h>
//Include this to show the errors in the exiting functions:
// #include <math.h>
//Include these to test the fixes:
#include "asinh.c"
#include "asinhf.c"
#include "asinhl.c"
#include "atanh.c"
#include "atanhf.c"
#include "atanhl.c"
void printErr (char *sMsg, int ndps, double x, double y) {
char xStr[50], yStr[50];
sprintf(xStr, "%.*e", ndps, x);
sprintf(yStr, "%.*e", ndps, y);
if (strcmp(xStr,yStr) != 0)
printf("%s: %.16e doesn't match %.16e\n",sMsg,x,y);
}
void test_asinh() {
double x;
x = 0.0; //use a variable to prevent constant folding (which uses a better
version of asinh).
printErr("asinh", 8, asinh( x), x);
printErr("asinh", 8, asinh(-x), -x);
for (x = -700.0; x <= 700.0; x += 10) {
printErr("asinh(sinh(x))", 12, asinh(sinh(x)),x);
}
for (x = -1; x <= 1; x += 0.01) {
printErr("asinh(sinh(x))", 12, asinh(sinh(x)),x);
}
for (x = -1e-300; x <= 1e-300; x += 1e-302) {
printErr("asinh(sinh(x))", 12, asinh(sinh(x)),x);
}
}
void test_asinhf() {
float xf;
xf = 0.0;
printErr("asinhf", 8, (double)asinhf( xf), xf);
printErr("asinhf", 8, (double)asinhf(-xf), -xf);
for (xf = -80.0; xf <= 80.0; xf += 1) {
printErr("asinhf(sinhf(xf))", 9, (double)asinhf(sinhf(xf)), (double)xf);
}
for (xf = -1; xf <= 1; xf += 0.01) {
//strangely float seems less accurate in the range +/- [0.5 - 1] (neither
sinh nor asinh).
printErr("asinhf(sinhf(xf))", 5, (double)asinhf(sinhf(xf)), (double)xf);
}
for (xf = -1e-30; xf <= 1e-30; xf += 1e-32) {
printErr("asinhf(sinhf(xf))", 9, (double)asinhf(sinhf(xf)), (double)xf);
}
}
void test_asinhl() {
long double xl;
xl = 0.0L;
printErr("asinhl", 8, (double)asinhl( xl), xl);
printErr("asinhl", 8, (double)asinhl(-xl), -xl);
for (xl = -11000.0L; xl <= 11000.0L; xl += 100.0L) {
//sinh(-10000L) is too big for a double, but asinh(sinh(-10000L)) is fine.
printErr("asinhl(sinhl(xl))", 15, (double)asinhl(sinhl(xl)), (double)xl);
}
for (xl = -1; xl <= 1; xl += 0.01) {
printErr("asinhl(sinhl(xl))", 15, (double)asinhl(sinhl(xl)), (double)xl);
}
for (xl = -1e-4900L; xl <= 1e-4900L; xl += 1e-4902L) {
printErr("asinhl(sinhl(xl))", 15, (double)asinhl(sinhl(xl)), (double)xl);
}
}
void test_atanh() {
double x;
x = 0.0;
printErr("atanh", 8, atanh( x), x);
printErr("atanh", 8, atanh(-x), -x);
/* only test tanh in range within [-pi, pi] to avoid infinities */
for (x = -1; x <= 1; x += 0.01) {
printErr("atanh(tanh(x))", 12, atanh(tanh(x)),x);
}
}
void test_atanhf() {
float xf;
xf = 0.0;
printErr("atanhf", 8, (double)atanhf( xf), xf);
printErr("atanhf", 8, (double)atanhf(-xf), -xf);
for (xf = -1; xf <= 1; xf += 0.01) {
printErr("atanhf(tanhf(xf))", 5, (double)atanhf(tanhf(xf)), (double)xf);
}
}
void test_atanhl() {
long double xl;
xl = 0.0L;
printErr("atanhl", 8, (double)atanhl( xl), xl);
printErr("atanhl", 8, (double)atanhl(-xl), -xl);
for (xl = -1; xl <= 1; xl += 0.01) {
printErr("atanhl(tanhl(xl))", 14, (double)atanhl(tanhl(xl)), (double)xl);
}
}
void main () {
test_asinh();
test_asinhf();
test_asinhl();
test_atanh();
test_atanhf();
test_atanhl();
printf("Completed tests.\n");
}
_______________________________________________
Mingw-w64-public mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/mingw-w64-public
