From: Paul Zimmermann <Paul.Zimmermann@inria.fr>
To: Szabolcs Nagy <nsz@port70.net>
Cc: musl@lists.openwall.com
Subject: [musl] Re: [PATCH 2/2] math: fix ld80 powl(x,huge) and powl(LDBL_MAX,small)
Date: Sat, 19 Aug 2023 08:14:28 +0200 [thread overview]
Message-ID: <p9u0v8db8sor.fsf@coriandre.loria.fr> (raw)
In-Reply-To: <20230818211738.GF3448312@port70.net> (message from Szabolcs Nagy on Fri, 18 Aug 2023 23:17:38 +0200)
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
>
>
prev parent reply other threads:[~2023-08-19 6:14 UTC|newest]
Thread overview: 6+ messages / expand[flat|nested] mbox.gz Atom feed top
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 ` Paul Zimmermann [this message]
Reply instructions:
You may reply publicly to this message via plain-text email
using any one of the following methods:
* Save the following mbox file, import it into your mail client,
and reply-to-all from there: mbox
Avoid top-posting and favor interleaved quoting:
https://en.wikipedia.org/wiki/Posting_style#Interleaved_style
* Reply using the --to, --cc, and --in-reply-to
switches of git-send-email(1):
git send-email \
--in-reply-to=p9u0v8db8sor.fsf@coriandre.loria.fr \
--to=paul.zimmermann@inria.fr \
--cc=musl@lists.openwall.com \
--cc=nsz@port70.net \
/path/to/YOUR_REPLY
https://kernel.org/pub/software/scm/git/docs/git-send-email.html
* If your mail client supports setting the In-Reply-To header
via mailto: links, try the mailto: link
Be sure your reply has a Subject: header at the top and a blank line
before the message body.
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).