From 8ec5b01346114da38e806ca1867da688d3a360e2 Mon Sep 17 00:00:00 2001
From: Joseph Myers <joseph@codesourcery.com>
Date: Sat, 29 Sep 2012 18:31:54 +0000
Subject: Fix sign of exact zero return from fma (bug 14638).

---
 sysdeps/ieee754/dbl-64/s_fma.c    | 5 +++++
 sysdeps/ieee754/dbl-64/s_fmaf.c   | 7 +++++++
 sysdeps/ieee754/ldbl-128/s_fma.c  | 8 +++++++-
 sysdeps/ieee754/ldbl-128/s_fmal.c | 5 +++++
 sysdeps/ieee754/ldbl-96/s_fma.c   | 6 +++++-
 sysdeps/ieee754/ldbl-96/s_fmal.c  | 5 +++++
 6 files changed, 34 insertions(+), 2 deletions(-)

(limited to 'sysdeps')

diff --git a/sysdeps/ieee754/dbl-64/s_fma.c b/sysdeps/ieee754/dbl-64/s_fma.c
index ce3bd36fac..c9809fb180 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 e7a0650f0f..a4f12d9f76 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);
diff --git a/sysdeps/ieee754/ldbl-128/s_fma.c b/sysdeps/ieee754/ldbl-128/s_fma.c
index 355b60ebbb..b08ff29c04 100644
--- a/sysdeps/ieee754/ldbl-128/s_fma.c
+++ b/sysdeps/ieee754/ldbl-128/s_fma.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.
 
@@ -33,6 +33,12 @@ __fma (double x, double y, double z)
   fenv_t env;
   /* Multiplication is always exact.  */
   long double temp = (long double) x * (long 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 (double) temp + z;
+
   union ieee854_long_double u;
   feholdexcept (&env);
   fesetround (FE_TOWARDZERO);
diff --git a/sysdeps/ieee754/ldbl-128/s_fmal.c b/sysdeps/ieee754/ldbl-128/s_fmal.c
index 963bbd7345..df68ade435 100644
--- a/sysdeps/ieee754/ldbl-128/s_fmal.c
+++ b/sysdeps/ieee754/ldbl-128/s_fmal.c
@@ -129,6 +129,11 @@ __fmal (long double x, long double y, long 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 ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1)
   long double x1 = x * C;
diff --git a/sysdeps/ieee754/ldbl-96/s_fma.c b/sysdeps/ieee754/ldbl-96/s_fma.c
index 78c0b0db18..001d8063d4 100644
--- a/sysdeps/ieee754/ldbl-96/s_fma.c
+++ b/sysdeps/ieee754/ldbl-96/s_fma.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.
 
@@ -38,6 +38,10 @@ __fma (double x, double y, double z)
       return (x * y) + z;
     }
 
+  /* 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 ((1ULL << (LDBL_MANT_DIG + 1) / 2) + 1)
   long double x1 = (long double) x * C;
diff --git a/sysdeps/ieee754/ldbl-96/s_fmal.c b/sysdeps/ieee754/ldbl-96/s_fmal.c
index ca1e0905a7..c27b0bd852 100644
--- a/sysdeps/ieee754/ldbl-96/s_fmal.c
+++ b/sysdeps/ieee754/ldbl-96/s_fmal.c
@@ -129,6 +129,11 @@ __fmal (long double x, long double y, long 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 ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1)
   long double x1 = x * C;
-- 
cgit 1.4.1