diff options
Diffstat (limited to 'libquadmath/math/fmaq.c')
| -rw-r--r-- | libquadmath/math/fmaq.c | 36 |
1 files changed, 20 insertions, 16 deletions
diff --git a/libquadmath/math/fmaq.c b/libquadmath/math/fmaq.c index 126b0a2d26b..23e3188669e 100644 --- a/libquadmath/math/fmaq.c +++ b/libquadmath/math/fmaq.c @@ -1,5 +1,5 @@ /* Compute x * y + z as ternary operation. - Copyright (C) 2010 Free Software Foundation, Inc. + Copyright (C) 2010-2012 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jakub@redhat.com>, 2010. @@ -57,6 +57,11 @@ fmaq (__float128 x, __float128 y, __float128 z) && u.ieee.exponent != 0x7fff && v.ieee.exponent != 0x7fff) return (z + x) + y; + /* If z is zero and x are y are nonzero, compute the result + as x * y to avoid the wrong sign of a zero result if x * y + underflows to 0. */ + if (z == 0 && x != 0 && y != 0) + return x * y; /* If x or y or z is Inf/NaN, or if fma will certainly overflow, or if x * y is less than half of FLT128_DENORM_MIN, compute as x * y + z. */ @@ -136,6 +141,11 @@ fmaq (__float128 x, __float128 y, __float128 z) y = v.value; z = w.value; } + + /* Ensure correct sign of exact 0 + 0. */ + if (__builtin_expect ((x == 0 || y == 0) && z == 0, 0)) + return x * y + z; + /* Multiplication m1 + m2 = x * y using Dekker's algorithm. */ #define C ((1LL << (FLT128_MANT_DIG + 1) / 2) + 1) __float128 x1 = x * C; @@ -191,7 +201,7 @@ fmaq (__float128 x, __float128 y, __float128 z) #endif v.value = a1 + u.value; /* Ensure the addition is not scheduled after fetestexcept call. */ - asm volatile ("" : : "m" (v)); + asm volatile ("" : : "m" (v.value)); #ifdef USE_FENV_H int j = fetestexcept (FE_INEXACT) != 0; feupdateenv (&env); @@ -220,20 +230,14 @@ fmaq (__float128 x, __float128 y, __float128 z) { /* v.ieee.mant_low & 2 is LSB bit of the result before rounding, v.ieee.mant_low & 1 is the round bit and j is our sticky - bit. In round-to-nearest 001 rounds down like 00, - 011 rounds up, even though 01 rounds down (thus we need - to adjust), 101 rounds down like 10 and 111 rounds up - like 11. */ - if ((v.ieee.mant_low & 3) == 1) - { - v.value *= 0x1p-226Q; - if (v.ieee.negative) - return v.value - 0x1p-16494Q /* __FLT128_DENORM_MIN__ */; - else - return v.value + 0x1p-16494Q /* __FLT128_DENORM_MIN__ */; - } - else - return v.value * 0x1p-226Q; + bit. */ + w.value = 0.0Q; + w.ieee.mant_low = ((v.ieee.mant_low & 3) << 1) | j; + w.ieee.negative = v.ieee.negative; + v.ieee.mant_low &= ~3U; + v.value *= 0x1p-226L; + w.value *= 0x1p-2L; + return v.value + w.value; } v.ieee.mant_low |= j; return v.value * 0x1p-226Q; |