diff options
| -rw-r--r-- | src/math/scalbn.c | 10 | ||||
| -rw-r--r-- | src/math/scalbnf.c | 8 | ||||
| -rw-r--r-- | 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; } |