Fix sign of exact zero return from fma (bug 14638).

diff --git a/ChangeLog b/ChangeLog
index 1d38d55029d4c4d61935e1b325443f12ff515f15..9e8fca0b0f3d1396266542c9f71515c91a28bf75 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,21 @@
+2012-09-29  Joseph Myers  <joseph@codesourcery.com>
+
+       [BZ #14638]
+       * sysdeps/ieee754/dbl-64/s_fma.c (__fma): Use x * y + z for exact
+       0 + 0.
+       * sysdeps/ieee754/dbl-64/s_fmaf.c (__fmaf): Use original rounding
+       mode for addition resulting in exact zero.
+       * sysdeps/ieee754/ldbl-128/s_fma.c (__fma): Likewise.
+       * sysdeps/ieee754/ldbl-128/s_fmal.c (__fmal): Use x * y + z for
+       exact 0 + 0.
+       * sysdeps/ieee754/ldbl-96/s_fma.c (__fma): Likewise.
+       * sysdeps/ieee754/ldbl-96/s_fmal.c (__fmal): Likewise.
+       * math/libm-test.inc (fma_test): Add more tests.
+       (fma_test_towardzero): New function.
+       (fma_test_downward): Likewise.
+       (fma_test_upward): Likewise.
+       (main): Call the new functions.
+
 2012-09-28  David S. Miller  <davem@davemloft.net>
 
        * sysdeps/sparc/fpu/libm-test-ulps: Fix garbage in file.

diff --git a/NEWS b/NEWS
index 0567c9e7073a242be4abd3f1d55c81dbbc965535..e816c23069307493766d5eacaa345ba5692cf6be 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -15,7 +15,7 @@ Version 2.17
   14195, 14237, 14252, 14283, 14298, 14303, 14307, 14328, 14331, 14336,
   14337, 14347, 14349, 14376, 14459, 14476, 14505, 14510, 14516, 14518,
   14519, 14530, 14532, 14538, 14543, 14544, 14545, 14562, 14576, 14579,
-  14583, 14587, 14621.
+  14583, 14587, 14621, 14638.
 
 * Support for STT_GNU_IFUNC symbols added for s390 and s390x.
   Optimized versions of memcpy, memset, and memcmp added for System z10 and

diff --git a/math/libm-test.inc b/math/libm-test.inc
index e8398bd0eea04411e4fecfa8ff034d4b18c9b9f7..007eea1f30658f10d440bda8ea3697e44ef3fa9c 100644 (file)
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -4546,6 +4546,36 @@ fma_test (void)
   TEST_fff_f (fma, minus_infty, minus_infty, plus_infty, plus_infty);
   TEST_fff_f (fma, plus_infty, minus_infty, minus_infty, minus_infty);
 
+  TEST_fff_f (fma, plus_zero, plus_zero, plus_zero, plus_zero);
+  TEST_fff_f (fma, plus_zero, plus_zero, minus_zero, plus_zero);
+  TEST_fff_f (fma, plus_zero, minus_zero, plus_zero, plus_zero);
+  TEST_fff_f (fma, plus_zero, minus_zero, minus_zero, minus_zero);
+  TEST_fff_f (fma, minus_zero, plus_zero, plus_zero, plus_zero);
+  TEST_fff_f (fma, minus_zero, plus_zero, minus_zero, minus_zero);
+  TEST_fff_f (fma, minus_zero, minus_zero, plus_zero, plus_zero);
+  TEST_fff_f (fma, minus_zero, minus_zero, minus_zero, plus_zero);
+  TEST_fff_f (fma, 1.0, plus_zero, plus_zero, plus_zero);
+  TEST_fff_f (fma, 1.0, plus_zero, minus_zero, plus_zero);
+  TEST_fff_f (fma, 1.0, minus_zero, plus_zero, plus_zero);
+  TEST_fff_f (fma, 1.0, minus_zero, minus_zero, minus_zero);
+  TEST_fff_f (fma, -1.0, plus_zero, plus_zero, plus_zero);
+  TEST_fff_f (fma, -1.0, plus_zero, minus_zero, minus_zero);
+  TEST_fff_f (fma, -1.0, minus_zero, plus_zero, plus_zero);
+  TEST_fff_f (fma, -1.0, minus_zero, minus_zero, plus_zero);
+  TEST_fff_f (fma, plus_zero, 1.0, plus_zero, plus_zero);
+  TEST_fff_f (fma, plus_zero, 1.0, minus_zero, plus_zero);
+  TEST_fff_f (fma, plus_zero, -1.0, plus_zero, plus_zero);
+  TEST_fff_f (fma, plus_zero, -1.0, minus_zero, minus_zero);
+  TEST_fff_f (fma, minus_zero, 1.0, plus_zero, plus_zero);
+  TEST_fff_f (fma, minus_zero, 1.0, minus_zero, minus_zero);
+  TEST_fff_f (fma, minus_zero, -1.0, plus_zero, plus_zero);
+  TEST_fff_f (fma, minus_zero, -1.0, minus_zero, plus_zero);
+
+  TEST_fff_f (fma, 1.0, 1.0, -1.0, plus_zero);
+  TEST_fff_f (fma, 1.0, -1.0, 1.0, plus_zero);
+  TEST_fff_f (fma, -1.0, 1.0, 1.0, plus_zero);
+  TEST_fff_f (fma, -1.0, -1.0, -1.0, plus_zero);
+
 #if defined (TEST_FLOAT) && FLT_MANT_DIG == 24
   TEST_fff_f (fma, 0x1.7ff8p+13, 0x1.000002p+0, 0x1.ffffp-24, 0x1.7ff802p+13);
   TEST_fff_f (fma, 0x1.fffp+0, 0x1.00001p+0, -0x1.fffp+0, 0x1.fffp-20);
@@ -4607,6 +4637,147 @@ fma_test (void)
 }
 
 
+static void
+fma_test_towardzero (void)
+{
+  int save_round_mode;
+  START (fma_towardzero);
+
+  save_round_mode = fegetround ();
+
+  if (!fesetround (FE_TOWARDZERO))
+    {
+      TEST_fff_f (fma, plus_zero, plus_zero, plus_zero, plus_zero);
+      TEST_fff_f (fma, plus_zero, plus_zero, minus_zero, plus_zero);
+      TEST_fff_f (fma, plus_zero, minus_zero, plus_zero, plus_zero);
+      TEST_fff_f (fma, plus_zero, minus_zero, minus_zero, minus_zero);
+      TEST_fff_f (fma, minus_zero, plus_zero, plus_zero, plus_zero);
+      TEST_fff_f (fma, minus_zero, plus_zero, minus_zero, minus_zero);
+      TEST_fff_f (fma, minus_zero, minus_zero, plus_zero, plus_zero);
+      TEST_fff_f (fma, minus_zero, minus_zero, minus_zero, plus_zero);
+      TEST_fff_f (fma, 1.0, plus_zero, plus_zero, plus_zero);
+      TEST_fff_f (fma, 1.0, plus_zero, minus_zero, plus_zero);
+      TEST_fff_f (fma, 1.0, minus_zero, plus_zero, plus_zero);
+      TEST_fff_f (fma, 1.0, minus_zero, minus_zero, minus_zero);
+      TEST_fff_f (fma, -1.0, plus_zero, plus_zero, plus_zero);
+      TEST_fff_f (fma, -1.0, plus_zero, minus_zero, minus_zero);
+      TEST_fff_f (fma, -1.0, minus_zero, plus_zero, plus_zero);
+      TEST_fff_f (fma, -1.0, minus_zero, minus_zero, plus_zero);
+      TEST_fff_f (fma, plus_zero, 1.0, plus_zero, plus_zero);
+      TEST_fff_f (fma, plus_zero, 1.0, minus_zero, plus_zero);
+      TEST_fff_f (fma, plus_zero, -1.0, plus_zero, plus_zero);
+      TEST_fff_f (fma, plus_zero, -1.0, minus_zero, minus_zero);
+      TEST_fff_f (fma, minus_zero, 1.0, plus_zero, plus_zero);
+      TEST_fff_f (fma, minus_zero, 1.0, minus_zero, minus_zero);
+      TEST_fff_f (fma, minus_zero, -1.0, plus_zero, plus_zero);
+      TEST_fff_f (fma, minus_zero, -1.0, minus_zero, plus_zero);
+
+      TEST_fff_f (fma, 1.0, 1.0, -1.0, plus_zero);
+      TEST_fff_f (fma, 1.0, -1.0, 1.0, plus_zero);
+      TEST_fff_f (fma, -1.0, 1.0, 1.0, plus_zero);
+      TEST_fff_f (fma, -1.0, -1.0, -1.0, plus_zero);
+    }
+
+  fesetround (save_round_mode);
+
+  END (fma_towardzero);
+}
+
+
+static void
+fma_test_downward (void)
+{
+  int save_round_mode;
+  START (fma_downward);
+
+  save_round_mode = fegetround ();
+
+  if (!fesetround (FE_DOWNWARD))
+    {
+      TEST_fff_f (fma, plus_zero, plus_zero, plus_zero, plus_zero);
+      TEST_fff_f (fma, plus_zero, plus_zero, minus_zero, minus_zero);
+      TEST_fff_f (fma, plus_zero, minus_zero, plus_zero, minus_zero);
+      TEST_fff_f (fma, plus_zero, minus_zero, minus_zero, minus_zero);
+      TEST_fff_f (fma, minus_zero, plus_zero, plus_zero, minus_zero);
+      TEST_fff_f (fma, minus_zero, plus_zero, minus_zero, minus_zero);
+      TEST_fff_f (fma, minus_zero, minus_zero, plus_zero, plus_zero);
+      TEST_fff_f (fma, minus_zero, minus_zero, minus_zero, minus_zero);
+      TEST_fff_f (fma, 1.0, plus_zero, plus_zero, plus_zero);
+      TEST_fff_f (fma, 1.0, plus_zero, minus_zero, minus_zero);
+      TEST_fff_f (fma, 1.0, minus_zero, plus_zero, minus_zero);
+      TEST_fff_f (fma, 1.0, minus_zero, minus_zero, minus_zero);
+      TEST_fff_f (fma, -1.0, plus_zero, plus_zero, minus_zero);
+      TEST_fff_f (fma, -1.0, plus_zero, minus_zero, minus_zero);
+      TEST_fff_f (fma, -1.0, minus_zero, plus_zero, plus_zero);
+      TEST_fff_f (fma, -1.0, minus_zero, minus_zero, minus_zero);
+      TEST_fff_f (fma, plus_zero, 1.0, plus_zero, plus_zero);
+      TEST_fff_f (fma, plus_zero, 1.0, minus_zero, minus_zero);
+      TEST_fff_f (fma, plus_zero, -1.0, plus_zero, minus_zero);
+      TEST_fff_f (fma, plus_zero, -1.0, minus_zero, minus_zero);
+      TEST_fff_f (fma, minus_zero, 1.0, plus_zero, minus_zero);
+      TEST_fff_f (fma, minus_zero, 1.0, minus_zero, minus_zero);
+      TEST_fff_f (fma, minus_zero, -1.0, plus_zero, plus_zero);
+      TEST_fff_f (fma, minus_zero, -1.0, minus_zero, minus_zero);
+
+      TEST_fff_f (fma, 1.0, 1.0, -1.0, minus_zero);
+      TEST_fff_f (fma, 1.0, -1.0, 1.0, minus_zero);
+      TEST_fff_f (fma, -1.0, 1.0, 1.0, minus_zero);
+      TEST_fff_f (fma, -1.0, -1.0, -1.0, minus_zero);
+    }
+
+  fesetround (save_round_mode);
+
+  END (fma_downward);
+}
+
+
+static void
+fma_test_upward (void)
+{
+  int save_round_mode;
+  START (fma_upward);
+
+  save_round_mode = fegetround ();
+
+  if (!fesetround (FE_UPWARD))
+    {
+      TEST_fff_f (fma, plus_zero, plus_zero, plus_zero, plus_zero);
+      TEST_fff_f (fma, plus_zero, plus_zero, minus_zero, plus_zero);
+      TEST_fff_f (fma, plus_zero, minus_zero, plus_zero, plus_zero);
+      TEST_fff_f (fma, plus_zero, minus_zero, minus_zero, minus_zero);
+      TEST_fff_f (fma, minus_zero, plus_zero, plus_zero, plus_zero);
+      TEST_fff_f (fma, minus_zero, plus_zero, minus_zero, minus_zero);
+      TEST_fff_f (fma, minus_zero, minus_zero, plus_zero, plus_zero);
+      TEST_fff_f (fma, minus_zero, minus_zero, minus_zero, plus_zero);
+      TEST_fff_f (fma, 1.0, plus_zero, plus_zero, plus_zero);
+      TEST_fff_f (fma, 1.0, plus_zero, minus_zero, plus_zero);
+      TEST_fff_f (fma, 1.0, minus_zero, plus_zero, plus_zero);
+      TEST_fff_f (fma, 1.0, minus_zero, minus_zero, minus_zero);
+      TEST_fff_f (fma, -1.0, plus_zero, plus_zero, plus_zero);
+      TEST_fff_f (fma, -1.0, plus_zero, minus_zero, minus_zero);
+      TEST_fff_f (fma, -1.0, minus_zero, plus_zero, plus_zero);
+      TEST_fff_f (fma, -1.0, minus_zero, minus_zero, plus_zero);
+      TEST_fff_f (fma, plus_zero, 1.0, plus_zero, plus_zero);
+      TEST_fff_f (fma, plus_zero, 1.0, minus_zero, plus_zero);
+      TEST_fff_f (fma, plus_zero, -1.0, plus_zero, plus_zero);
+      TEST_fff_f (fma, plus_zero, -1.0, minus_zero, minus_zero);
+      TEST_fff_f (fma, minus_zero, 1.0, plus_zero, plus_zero);
+      TEST_fff_f (fma, minus_zero, 1.0, minus_zero, minus_zero);
+      TEST_fff_f (fma, minus_zero, -1.0, plus_zero, plus_zero);
+      TEST_fff_f (fma, minus_zero, -1.0, minus_zero, plus_zero);
+
+      TEST_fff_f (fma, 1.0, 1.0, -1.0, plus_zero);
+      TEST_fff_f (fma, 1.0, -1.0, 1.0, plus_zero);
+      TEST_fff_f (fma, -1.0, 1.0, 1.0, plus_zero);
+      TEST_fff_f (fma, -1.0, -1.0, -1.0, plus_zero);
+    }
+
+  fesetround (save_round_mode);
+
+  END (fma_upward);
+}
+
+
 static void
 fmax_test (void)
 {
@@ -9539,6 +9710,9 @@ main (int argc, char **argv)
 
   /* Multiply and add:  */
   fma_test ();
+  fma_test_towardzero ();
+  fma_test_downward ();
+  fma_test_upward ();
 
   /* Complex functions:  */
   cabs_test ();

diff --git a/sysdeps/ieee754/dbl-64/s_fma.c b/sysdeps/ieee754/dbl-64/s_fma.c
index ce3bd36fac4eebcf956f469661d5c81397e48a83..c9809fb180685a941044f4907da4a275f69e3c6d 100644 (file)
--- 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 e7a0650f0f3036beb3f6d164bb0ea7fe1222bb58..a4f12d9f76fa73d4f4cd926be414df123b4b74f5 100644 (file)
--- 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 355b60ebbbedb0543113390a5abbf84cf9710257..b08ff29c0408cd259c1573ea94af5d1ab2c22e5e 100644 (file)
--- 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 963bbd73454f8e34f92431183726c7b8ff8957cd..df68ade435f2707856c0b7e0cca8ffb8e6e4dd7d 100644 (file)
--- 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 78c0b0db18174767d02b3d11284add1ec960496a..001d8063d4a9c6f5581d4d8d685682afa7d7443f 100644 (file)
--- 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 ca1e0905a7997e02004732cc8166fbaaabb44e6b..c27b0bd852c92ee7de6212ddbfa4551724373c03 100644 (file)
--- 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;