diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64')
-rw-r--r-- | sysdeps/ieee754/dbl-64/s_fma.c | 5 | ||||
-rw-r--r-- | sysdeps/ieee754/dbl-64/s_fmaf.c | 7 |
2 files changed, 12 insertions, 0 deletions
diff --git a/sysdeps/ieee754/dbl-64/s_fma.c b/sysdeps/ieee754/dbl-64/s_fma.c index ce3bd36..c9809fb 100644 --- a/sysdeps/ieee754/dbl-64/s_fma.c +++ b/sysdeps/ieee754/dbl-64/s_fma.c @@ -128,6 +128,11 @@ __fma (double x, double y, double z) y = v.d; z = w.d; } + + /* 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 ((1 << (DBL_MANT_DIG + 1) / 2) + 1) double x1 = x * C; diff --git a/sysdeps/ieee754/dbl-64/s_fmaf.c b/sysdeps/ieee754/dbl-64/s_fmaf.c index e7a0650..a4f12d9 100644 --- a/sysdeps/ieee754/dbl-64/s_fmaf.c +++ b/sysdeps/ieee754/dbl-64/s_fmaf.c @@ -32,8 +32,15 @@ float __fmaf (float x, float y, float z) { fenv_t env; + /* Multiplication is always exact. */ double temp = (double) x * (double) y; + + /* Ensure correct sign of an exact zero result by performing the + addition in the original rounding mode in that case. */ + if (temp == -z) + return (float) temp + z; + union ieee754_double u; libc_feholdexcept_setround (&env, FE_TOWARDZERO); |