RepoMirrors/musl
1
0
Fork 0
mirror of git://git.musl-libc.org/musl synced 2025-03-04 10:47:32 +00:00
Code Issues Packages Projects Releases Wiki Activity

math: fix ld80 powl(x,huge) and powl(LDBL_MAX,small)

Browse Source 
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).
This commit is contained in:
Szabolcs Nagy 2023-08-18 23:17:38 +02:00 committed by Rich Felker
parent 6d10102709
commit 39e43f0881
1 changed files with 21 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

34
src/math/powl.c
View File

@ -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)
/* 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 (y <= -LDBL_MAX) {
if (x > 1.0 || x < -1.0)
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;