diff options
Diffstat (limited to 'sysdeps/ieee754/ldbl-128/s_fmal.c')
| -rw-r--r-- | sysdeps/ieee754/ldbl-128/s_fmal.c | 45 |
1 files changed, 27 insertions, 18 deletions
diff --git a/sysdeps/ieee754/ldbl-128/s_fmal.c b/sysdeps/ieee754/ldbl-128/s_fmal.c index 6fa663a6c0..c9accad8a3 100644 --- a/sysdeps/ieee754/ldbl-128/s_fmal.c +++ b/sysdeps/ieee754/ldbl-128/s_fmal.c @@ -118,8 +118,17 @@ __fmal (long double x, long double y, long double z) { /* Similarly. If z exponent is very large and x and y exponents are - very small, it doesn't matter if we don't adjust it. */ - if (u.ieee.exponent > v.ieee.exponent) + very small, adjust them up to avoid spurious underflows, + rather than down. */ + if (u.ieee.exponent + v.ieee.exponent + <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG) + { + if (u.ieee.exponent > v.ieee.exponent) + u.ieee.exponent += 2 * LDBL_MANT_DIG + 2; + else + v.ieee.exponent += 2 * LDBL_MANT_DIG + 2; + } + else if (u.ieee.exponent > v.ieee.exponent) { if (u.ieee.exponent > LDBL_MANT_DIG) u.ieee.exponent -= LDBL_MANT_DIG; @@ -149,15 +158,15 @@ __fmal (long double x, long double y, long double z) <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG) */ { if (u.ieee.exponent > v.ieee.exponent) - u.ieee.exponent += 2 * LDBL_MANT_DIG; + u.ieee.exponent += 2 * LDBL_MANT_DIG + 2; else - v.ieee.exponent += 2 * LDBL_MANT_DIG; - if (w.ieee.exponent <= 4 * LDBL_MANT_DIG + 4) + v.ieee.exponent += 2 * LDBL_MANT_DIG + 2; + if (w.ieee.exponent <= 4 * LDBL_MANT_DIG + 6) { if (w.ieee.exponent) - w.ieee.exponent += 2 * LDBL_MANT_DIG; + w.ieee.exponent += 2 * LDBL_MANT_DIG + 2; else - w.d *= 0x1p226L; + w.d *= 0x1p228L; adjust = -1; } /* Otherwise x * y should just affect inexact @@ -242,19 +251,19 @@ __fmal (long double x, long double y, long double z) /* If a1 + u.d is exact, the only rounding happens during scaling down. */ if (j == 0) - return v.d * 0x1p-226L; + return v.d * 0x1p-228L; /* If result rounded to zero is not subnormal, no double rounding will occur. */ - if (v.ieee.exponent > 226) - return (a1 + u.d) * 0x1p-226L; - /* If v.d * 0x1p-226L with round to zero is a subnormal above - or equal to LDBL_MIN / 2, then v.d * 0x1p-226L shifts mantissa + if (v.ieee.exponent > 228) + return (a1 + u.d) * 0x1p-228L; + /* If v.d * 0x1p-228L with round to zero is a subnormal above + or equal to LDBL_MIN / 2, then v.d * 0x1p-228L shifts mantissa down just by 1 bit, which means v.ieee.mantissa3 |= j would change the round bit, not sticky or guard bit. - v.d * 0x1p-226L never normalizes by shifting up, + v.d * 0x1p-228L never normalizes by shifting up, so round bit plus sticky bit should be already enough for proper rounding. */ - if (v.ieee.exponent == 226) + if (v.ieee.exponent == 228) { /* If the exponent would be in the normal range when rounding to normal precision with unbounded exponent @@ -264,8 +273,8 @@ __fmal (long double x, long double y, long double z) if (TININESS_AFTER_ROUNDING) { w.d = a1 + u.d; - if (w.ieee.exponent == 227) - return w.d * 0x1p-226L; + if (w.ieee.exponent == 229) + return w.d * 0x1p-228L; } /* v.ieee.mantissa3 & 2 is LSB bit of the result before rounding, v.ieee.mantissa3 & 1 is the round bit and j is our sticky @@ -274,12 +283,12 @@ __fmal (long double x, long double y, long double z) w.ieee.mantissa3 = ((v.ieee.mantissa3 & 3) << 1) | j; w.ieee.negative = v.ieee.negative; v.ieee.mantissa3 &= ~3U; - v.d *= 0x1p-226L; + v.d *= 0x1p-228L; w.d *= 0x1p-2L; return v.d + w.d; } v.ieee.mantissa3 |= j; - return v.d * 0x1p-226L; + return v.d * 0x1p-228L; } } weak_alias (__fmal, fmal) |