* [musl] musl 1.2.4 @ 2023-08-16 14:39 Paul Zimmermann 2023-08-17 15:57 ` Szabolcs Nagy ` (2 more replies) 0 siblings, 3 replies; 6+ messages in thread From: Paul Zimmermann @ 2023-08-16 14:39 UTC (permalink / raw) To: musl Hi, while updating our comparison of mathematical functions with Musl 1.2.4, I noticed the following two issues with acoshl and exp10l in double extended precision: zimmerma@coriandre:~/svn/tbd/20/src/binary80$ VERBOSE=-v ./doit.musl acosh 1000 Checking acosh with musl-1.2.4 Using seed 2804715 Using 6 threads NEW acosh 0 -1 -0x6.e2368c0ed74e5698p+16l [-nan] [inf] inf inf libm gives -0x4.b4d6a621e8e631f8p+0l mpfr gives nanl zimmerma@coriandre:~/svn/tbd/20/src/binary80$ VERBOSE=-v ./doit.musl exp10 1000 Checking exp10 with musl-1.2.4 Using seed 2807610 Using 6 threads exp10 0 -1 0x2.68826a13ef3fde64p+16376l [-nan] [inf] inf inf libm gives nanl mpfr gives infl These issues were already present in 1.2.2 at least, but only detected recently by our program. Best regards, Paul ^ permalink raw reply [flat|nested] 6+ messages in thread
* Re: [musl] musl 1.2.4 2023-08-16 14:39 [musl] musl 1.2.4 Paul Zimmermann @ 2023-08-17 15:57 ` Szabolcs Nagy 2023-08-18 21:16 ` [musl] [PATCH 1/2] math: fix ld80 acoshl(x) for x < 0 Szabolcs Nagy 2023-08-18 21:17 ` [musl] [PATCH 2/2] math: fix ld80 powl(x,huge) and powl(LDBL_MAX,small) Szabolcs Nagy 2 siblings, 0 replies; 6+ messages in thread From: Szabolcs Nagy @ 2023-08-17 15:57 UTC (permalink / raw) To: Paul Zimmermann; +Cc: musl * Paul Zimmermann <Paul.Zimmermann@inria.fr> [2023-08-16 16:39:25 +0200]: > while updating our comparison of mathematical functions with Musl 1.2.4, > I noticed the following two issues with acoshl and exp10l in double extended > precision: > > zimmerma@coriandre:~/svn/tbd/20/src/binary80$ VERBOSE=-v ./doit.musl acosh 1000 > Checking acosh with musl-1.2.4 > Using seed 2804715 > Using 6 threads > NEW acosh 0 -1 -0x6.e2368c0ed74e5698p+16l [-nan] [inf] inf inf > libm gives -0x4.b4d6a621e8e631f8p+0l > mpfr gives nanl yes it seems acoshl(x) does not handle x<0 (and can go wrong for x in [-2,-0x1p32]) > > zimmerma@coriandre:~/svn/tbd/20/src/binary80$ VERBOSE=-v ./doit.musl exp10 1000 > Checking exp10 with musl-1.2.4 > Using seed 2807610 > Using 6 threads > exp10 0 -1 0x2.68826a13ef3fde64p+16376l [-nan] [inf] inf inf > libm gives nanl > mpfr gives infl this seems to be a bug in powl: the reducl function can overflow. thanks for the reports. ^ permalink raw reply [flat|nested] 6+ messages in thread
* [musl] [PATCH 1/2] math: fix ld80 acoshl(x) for x < 0 2023-08-16 14:39 [musl] musl 1.2.4 Paul Zimmermann 2023-08-17 15:57 ` Szabolcs Nagy @ 2023-08-18 21:16 ` Szabolcs Nagy 2023-08-19 6:10 ` [musl] " Paul Zimmermann 2023-08-18 21:17 ` [musl] [PATCH 2/2] math: fix ld80 powl(x,huge) and powl(LDBL_MAX,small) Szabolcs Nagy 2 siblings, 1 reply; 6+ messages in thread From: Szabolcs Nagy @ 2023-08-18 21:16 UTC (permalink / raw) To: Paul Zimmermann; +Cc: musl acosh(x) is nan for x < 1, but x < 0 cases were not handled specially and acoshl gave wrong result for some -0x1p32 < x < -2 values, e.g.: acoshl(-0x1p20) returned -inf, acoshl(-0x1.4p20) returned -0x1.db365758403aa9acp+0L, fixed by checking the sign bit and handling it specially. reported by Paul Zimmermann. --- src/math/acoshl.c | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/math/acoshl.c b/src/math/acoshl.c index 8d4b43f6..943cec17 100644 --- a/src/math/acoshl.c +++ b/src/math/acoshl.c @@ -10,14 +10,18 @@ long double acoshl(long double x) long double acoshl(long double x) { union ldshape u = {x}; - int e = u.i.se & 0x7fff; + int e = u.i.se; if (e < 0x3fff + 1) - /* |x| < 2, invalid if x < 1 or nan */ + /* 0 <= x < 2, invalid if x < 1 */ return log1pl(x-1 + sqrtl((x-1)*(x-1)+2*(x-1))); if (e < 0x3fff + 32) - /* |x| < 0x1p32 */ + /* 2 <= x < 0x1p32 */ return logl(2*x - 1/(x+sqrtl(x*x-1))); + if (e & 0x8000) + /* x < 0 or x = -0, invalid */ + return (x - x) / (x - x); + /* 0x1p32 <= x or nan */ return logl(x) + 0.693147180559945309417232121458176568L; } #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 -- 2.41.0 ^ permalink raw reply [flat|nested] 6+ messages in thread
* [musl] Re: [PATCH 1/2] math: fix ld80 acoshl(x) for x < 0 2023-08-18 21:16 ` [musl] [PATCH 1/2] math: fix ld80 acoshl(x) for x < 0 Szabolcs Nagy @ 2023-08-19 6:10 ` Paul Zimmermann 0 siblings, 0 replies; 6+ messages in thread From: Paul Zimmermann @ 2023-08-19 6:10 UTC (permalink / raw) To: Szabolcs Nagy; +Cc: musl Dear Szabolcs, ok for me, I now get a largest error of 2.99 ulp: NEW acosh 0 -1 0x1.1ecdb5b8f0c5d79p+0l [3] [2.99] 2.98085 2.980840623325726 Thank you for fixing that rapidly! Paul > Date: Fri, 18 Aug 2023 23:16:00 +0200 > From: Szabolcs Nagy <nsz@port70.net> > Cc: musl@lists.openwall.com > > acosh(x) is nan for x < 1, but x < 0 cases were not handled specially > and acoshl gave wrong result for some -0x1p32 < x < -2 values, e.g.: > > acoshl(-0x1p20) returned -inf, > acoshl(-0x1.4p20) returned -0x1.db365758403aa9acp+0L, > > fixed by checking the sign bit and handling it specially. > > reported by Paul Zimmermann. > --- > src/math/acoshl.c | 10 +++++++--- > 1 file changed, 7 insertions(+), 3 deletions(-) > > diff --git a/src/math/acoshl.c b/src/math/acoshl.c > index 8d4b43f6..943cec17 100644 > --- a/src/math/acoshl.c > +++ b/src/math/acoshl.c > @@ -10,14 +10,18 @@ long double acoshl(long double x) > long double acoshl(long double x) > { > union ldshape u = {x}; > - int e = u.i.se & 0x7fff; > + int e = u.i.se; > > if (e < 0x3fff + 1) > - /* |x| < 2, invalid if x < 1 or nan */ > + /* 0 <= x < 2, invalid if x < 1 */ > return log1pl(x-1 + sqrtl((x-1)*(x-1)+2*(x-1))); > if (e < 0x3fff + 32) > - /* |x| < 0x1p32 */ > + /* 2 <= x < 0x1p32 */ > return logl(2*x - 1/(x+sqrtl(x*x-1))); > + if (e & 0x8000) > + /* x < 0 or x = -0, invalid */ > + return (x - x) / (x - x); > + /* 0x1p32 <= x or nan */ > return logl(x) + 0.693147180559945309417232121458176568L; > } > #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 > -- > 2.41.0 > > ^ permalink raw reply [flat|nested] 6+ messages in thread
* [musl] [PATCH 2/2] math: fix ld80 powl(x,huge) and powl(LDBL_MAX,small) 2023-08-16 14:39 [musl] musl 1.2.4 Paul Zimmermann 2023-08-17 15:57 ` Szabolcs Nagy 2023-08-18 21:16 ` [musl] [PATCH 1/2] math: fix ld80 acoshl(x) for x < 0 Szabolcs Nagy @ 2023-08-18 21:17 ` Szabolcs Nagy 2023-08-19 6:14 ` [musl] " Paul Zimmermann 2 siblings, 1 reply; 6+ messages in thread From: Szabolcs Nagy @ 2023-08-18 21:17 UTC (permalink / raw) To: Paul Zimmermann; +Cc: musl powl used >= LDBL_MAX as infinity check, but LDBL_MAX is finite, so this can cause wrong results e.g. powl(LDBL_MAX, 0.5) returned inf or powl(2, LDBL_MAX) returned inf without raising overflow. huge y values (close to LDBL_MAX) could cause intermediate results to overflow (computing y * log2(x) with more than long double precision) and e.g. powl(0.5, 0x1p16380L) or powl(10, 0x1p16380L) returned nan. this is fixed by handling huge y early since that always overflows or underflows. reported by Paul Zimmermann against expl10 (which uses powl). --- src/math/powl.c | 34 +++++++++++++++++++++------------- 1 file changed, 21 insertions(+), 13 deletions(-) diff --git a/src/math/powl.c b/src/math/powl.c index 5b6da07b..6f64ea71 100644 --- a/src/math/powl.c +++ b/src/math/powl.c @@ -212,25 +212,33 @@ long double powl(long double x, long double y) } if (x == 1.0) return 1.0; /* 1**y = 1, even if y is nan */ - if (x == -1.0 && !isfinite(y)) - return 1.0; /* -1**inf = 1 */ if (y == 0.0) return 1.0; /* x**0 = 1, even if x is nan */ if (y == 1.0) return x; - if (y >= LDBL_MAX) { - if (x > 1.0 || x < -1.0) - return INFINITY; - if (x != 0.0) - return 0.0; - } - if (y <= -LDBL_MAX) { - if (x > 1.0 || x < -1.0) + /* if y*log2(x) < log2(LDBL_TRUE_MIN)-1 then x^y uflows to 0 + if y*log2(x) > -log2(LDBL_TRUE_MIN)+1 > LDBL_MAX_EXP then x^y oflows + if |x|!=1 then |log2(x)| > |log(x)| > LDBL_EPSILON/2 so + x^y oflows/uflows if |y|*LDBL_EPSILON/2 > -log2(LDBL_TRUE_MIN)+1 */ + if (fabsl(y) > 2*(-LDBL_MIN_EXP+LDBL_MANT_DIG+1)/LDBL_EPSILON) { + /* y is not an odd int */ + if (x == -1.0) + return 1.0; + if (y == INFINITY) { + if (x > 1.0 || x < -1.0) + return INFINITY; return 0.0; - if (x != 0.0 || y == -INFINITY) + } + if (y == -INFINITY) { + if (x > 1.0 || x < -1.0) + return 0.0; return INFINITY; + } + if ((x > 1.0 || x < -1.0) == (y > 0)) + return huge * huge; + return twom10000 * twom10000; } - if (x >= LDBL_MAX) { + if (x == INFINITY) { if (y > 0.0) return INFINITY; return 0.0; @@ -253,7 +261,7 @@ long double powl(long double x, long double y) yoddint = 1; } - if (x <= -LDBL_MAX) { + if (x == -INFINITY) { if (y > 0.0) { if (yoddint) return -INFINITY; -- 2.41.0 ^ permalink raw reply [flat|nested] 6+ messages in thread
* [musl] Re: [PATCH 2/2] math: fix ld80 powl(x,huge) and powl(LDBL_MAX,small) 2023-08-18 21:17 ` [musl] [PATCH 2/2] math: fix ld80 powl(x,huge) and powl(LDBL_MAX,small) Szabolcs Nagy @ 2023-08-19 6:14 ` Paul Zimmermann 0 siblings, 0 replies; 6+ messages in thread From: Paul Zimmermann @ 2023-08-19 6:14 UTC (permalink / raw) To: Szabolcs Nagy; +Cc: musl Dear Szabolcs, looks good to me, I now get a largest error of 40.1 ulps for exp10l: NEW exp10 0 0 0xd.4135bbaf96d9716p+8l [40] [40.1] 40.0525 40.05240557189884 Paul > Date: Fri, 18 Aug 2023 23:17:38 +0200 > From: Szabolcs Nagy <nsz@port70.net> > Cc: musl@lists.openwall.com > > powl used >= LDBL_MAX as infinity check, but LDBL_MAX is finite, so > this can cause wrong results e.g. powl(LDBL_MAX, 0.5) returned inf > or powl(2, LDBL_MAX) returned inf without raising overflow. > > huge y values (close to LDBL_MAX) could cause intermediate results to > overflow (computing y * log2(x) with more than long double precision) > and e.g. powl(0.5, 0x1p16380L) or powl(10, 0x1p16380L) returned nan. > this is fixed by handling huge y early since that always overflows or > underflows. > > reported by Paul Zimmermann against expl10 (which uses powl). > --- > src/math/powl.c | 34 +++++++++++++++++++++------------- > 1 file changed, 21 insertions(+), 13 deletions(-) > > diff --git a/src/math/powl.c b/src/math/powl.c > index 5b6da07b..6f64ea71 100644 > --- a/src/math/powl.c > +++ b/src/math/powl.c > @@ -212,25 +212,33 @@ long double powl(long double x, long double y) > } > if (x == 1.0) > return 1.0; /* 1**y = 1, even if y is nan */ > - if (x == -1.0 && !isfinite(y)) > - return 1.0; /* -1**inf = 1 */ > if (y == 0.0) > return 1.0; /* x**0 = 1, even if x is nan */ > if (y == 1.0) > return x; > - if (y >= LDBL_MAX) { > - if (x > 1.0 || x < -1.0) > - return INFINITY; > - if (x != 0.0) > - return 0.0; > - } > - if (y <= -LDBL_MAX) { > - if (x > 1.0 || x < -1.0) > + /* if y*log2(x) < log2(LDBL_TRUE_MIN)-1 then x^y uflows to 0 > + if y*log2(x) > -log2(LDBL_TRUE_MIN)+1 > LDBL_MAX_EXP then x^y oflows > + if |x|!=1 then |log2(x)| > |log(x)| > LDBL_EPSILON/2 so > + x^y oflows/uflows if |y|*LDBL_EPSILON/2 > -log2(LDBL_TRUE_MIN)+1 */ > + if (fabsl(y) > 2*(-LDBL_MIN_EXP+LDBL_MANT_DIG+1)/LDBL_EPSILON) { > + /* y is not an odd int */ > + if (x == -1.0) > + return 1.0; > + if (y == INFINITY) { > + if (x > 1.0 || x < -1.0) > + return INFINITY; > return 0.0; > - if (x != 0.0 || y == -INFINITY) > + } > + if (y == -INFINITY) { > + if (x > 1.0 || x < -1.0) > + return 0.0; > return INFINITY; > + } > + if ((x > 1.0 || x < -1.0) == (y > 0)) > + return huge * huge; > + return twom10000 * twom10000; > } > - if (x >= LDBL_MAX) { > + if (x == INFINITY) { > if (y > 0.0) > return INFINITY; > return 0.0; > @@ -253,7 +261,7 @@ long double powl(long double x, long double y) > yoddint = 1; > } > > - if (x <= -LDBL_MAX) { > + if (x == -INFINITY) { > if (y > 0.0) { > if (yoddint) > return -INFINITY; > -- > 2.41.0 > > ^ permalink raw reply [flat|nested] 6+ messages in thread
end of thread, other threads:[~2023-08-19 6:14 UTC | newest] Thread overview: 6+ messages (download: mbox.gz / follow: Atom feed) -- links below jump to the message on this page -- 2023-08-16 14:39 [musl] musl 1.2.4 Paul Zimmermann 2023-08-17 15:57 ` Szabolcs Nagy 2023-08-18 21:16 ` [musl] [PATCH 1/2] math: fix ld80 acoshl(x) for x < 0 Szabolcs Nagy 2023-08-19 6:10 ` [musl] " Paul Zimmermann 2023-08-18 21:17 ` [musl] [PATCH 2/2] math: fix ld80 powl(x,huge) and powl(LDBL_MAX,small) Szabolcs Nagy 2023-08-19 6:14 ` [musl] " Paul Zimmermann
Code repositories for project(s) associated with this public inbox https://git.vuxu.org/mirror/musl/ This is a public inbox, see mirroring instructions for how to clone and mirror all data and code used for this inbox; as well as URLs for NNTP newsgroup(s).