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.)

diff --git a/src/math/scalbn.c b/src/math/scalbn.c
index 530e07c79f051aab246f39b31bfa1571699a8f57..182f561068fda7bc8321dfc7991df8b6d58e1b56 100644 (file)
--- 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 0b62c3c71f85d4b72a21a9b39752fb77ce0caf10..a5ad208b69929f24336fae40e6af37101c99a72f 100644 (file)
--- 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 08a4c58754d8a237f4441f10b1cce610ef885e43..db44dab06482748a06808903c9b011740feabf92 100644 (file)
--- 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;
                }