* [PATCH] fix scalbn when result is in the subnormal range
@ 2017-04-03 0:38 Szabolcs Nagy
2017-04-21 21:37 ` Rich Felker
0 siblings, 1 reply; 2+ messages in thread
From: Szabolcs Nagy @ 2017-04-03 0:38 UTC (permalink / raw)
To: musl
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.)
---
src/math/scalbn.c | 10 ++++++----
src/math/scalbnf.c | 8 ++++----
src/math/scalbnl.c | 8 ++++----
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;
}
--
2.12.2
^ permalink raw reply [flat|nested] 2+ messages in thread
* Re: [PATCH] fix scalbn when result is in the subnormal range
2017-04-03 0:38 [PATCH] fix scalbn when result is in the subnormal range Szabolcs Nagy
@ 2017-04-21 21:37 ` Rich Felker
0 siblings, 0 replies; 2+ messages in thread
From: Rich Felker @ 2017-04-21 21:37 UTC (permalink / raw)
To: musl
On Mon, Apr 03, 2017 at 02:38:13AM +0200, Szabolcs Nagy wrote:
> 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.)
OK, applying.
Rich
^ permalink raw reply [flat|nested] 2+ messages in thread
end of thread, other threads:[~2017-04-21 21:37 UTC | newest]
Thread overview: 2+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2017-04-03 0:38 [PATCH] fix scalbn when result is in the subnormal range Szabolcs Nagy
2017-04-21 21:37 ` Rich Felker
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).