[<prev] [next>] [<thread-prev] [day] [month] [year] [list]
Message-Id: <p9u0v8db8sor.fsf@coriandre.loria.fr>
Date: Sat, 19 Aug 2023 08:14:28 +0200
From: Paul Zimmermann <Paul.Zimmermann@...ia.fr>
To: Szabolcs Nagy <nsz@...t70.net>
Cc: musl@...ts.openwall.com
Subject: Re: [PATCH 2/2] math: fix ld80 powl(x,huge) and powl(LDBL_MAX,small)
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@...t70.net>
> Cc: musl@...ts.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
>
>
Powered by blists - more mailing lists
Confused about mailing lists and their use?
Read about mailing lists on Wikipedia
and check out these
guidelines on proper formatting of your messages.
