diff options
| author | Szabolcs Nagy <nsz@port70.net> | 2023-08-18 23:17:38 +0200 |
|---|---|---|
| committer | Rich Felker <dalias@aerifal.cx> | 2023-08-19 14:09:00 -0400 |
| commit | 39e43f0881b12d602ddcb8a9e9295456354b765b (patch) | |
| tree | 91a5db902649322e8700f0fdc4fc529418fbce1c /src/math | |
| parent | 6d10102709df4bc966d2846c1c45cd667e5048e5 (diff) | |
| download | musl-39e43f0881b12d602ddcb8a9e9295456354b765b.tar.gz musl-39e43f0881b12d602ddcb8a9e9295456354b765b.tar.xz musl-39e43f0881b12d602ddcb8a9e9295456354b765b.zip | |
math: fix ld80 powl(x,huge) and powl(LDBL_MAX,small)
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).
Diffstat (limited to 'src/math')
| -rw-r--r-- | src/math/powl.c | 34 |
1 files 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; |