summaryrefslogtreecommitdiff
path: root/src/math
diff options
context:
space:
mode:
authorSzabolcs Nagy <nsz@port70.net>2017-04-03 02:38:13 +0200
committerRich Felker <dalias@aerifal.cx>2017-04-21 17:31:00 -0400
commit8c44a060243f04283ca68dad199aab90336141db (patch)
tree5fd8fd700ac059c1ef2f2999c5783207c91bcf59 /src/math
parent2577b1bc16124d0690b9dd268a9f582f80bdcd67 (diff)
downloadmusl-8c44a060243f04283ca68dad199aab90336141db.tar.gz
musl-8c44a060243f04283ca68dad199aab90336141db.tar.bz2
musl-8c44a060243f04283ca68dad199aab90336141db.tar.xz
musl-8c44a060243f04283ca68dad199aab90336141db.zip
fix scalbn when result is in the subnormal range
in nearest rounding mode scalbn could introduce double rounding error when an intermediate value and the final result were both in the subnormal range e.g. scalbn(0x1.7ffffffffffffp-1, -1073) returned 0x1p-1073 instead of 0x1p-1074, because the intermediate computation got rounded to 0x1.8p-1023. with the fix an intermediate value can only be in the subnormal range if the final result is 0 which is correct even after double rounding. (there still can be two roundings so signals may be raised twice, but that's only observable with trapping exceptions which is not supported.)
Diffstat (limited to 'src/math')
-rw-r--r--src/math/scalbn.c10
-rw-r--r--src/math/scalbnf.c8
-rw-r--r--src/math/scalbnl.c8
3 files changed, 14 insertions, 12 deletions
diff --git a/src/math/scalbn.c b/src/math/scalbn.c
index 530e07c7..182f5610 100644
--- a/src/math/scalbn.c
+++ b/src/math/scalbn.c
@@ -16,11 +16,13 @@ double scalbn(double x, int n)
n = 1023;
}
} else if (n < -1022) {
- y *= 0x1p-1022;
- n += 1022;
+ /* make sure final n < -53 to avoid double
+ rounding in the subnormal range */
+ y *= 0x1p-1022 * 0x1p53;
+ n += 1022 - 53;
if (n < -1022) {
- y *= 0x1p-1022;
- n += 1022;
+ y *= 0x1p-1022 * 0x1p53;
+ n += 1022 - 53;
if (n < -1022)
n = -1022;
}
diff --git a/src/math/scalbnf.c b/src/math/scalbnf.c
index 0b62c3c7..a5ad208b 100644
--- a/src/math/scalbnf.c
+++ b/src/math/scalbnf.c
@@ -16,11 +16,11 @@ float scalbnf(float x, int n)
n = 127;
}
} else if (n < -126) {
- y *= 0x1p-126f;
- n += 126;
+ y *= 0x1p-126f * 0x1p24f;
+ n += 126 - 24;
if (n < -126) {
- y *= 0x1p-126f;
- n += 126;
+ y *= 0x1p-126f * 0x1p24f;
+ n += 126 - 24;
if (n < -126)
n = -126;
}
diff --git a/src/math/scalbnl.c b/src/math/scalbnl.c
index 08a4c587..db44dab0 100644
--- a/src/math/scalbnl.c
+++ b/src/math/scalbnl.c
@@ -20,11 +20,11 @@ long double scalbnl(long double x, int n)
n = 16383;
}
} else if (n < -16382) {
- x *= 0x1p-16382L;
- n += 16382;
+ x *= 0x1p-16382L * 0x1p113L;
+ n += 16382 - 113;
if (n < -16382) {
- x *= 0x1p-16382L;
- n += 16382;
+ x *= 0x1p-16382L * 0x1p113L;
+ n += 16382 - 113;
if (n < -16382)
n = -16382;
}