This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix sign of exact zero return from fma (bug 14638)
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: <libc-alpha at sourceware dot org>
- Cc: David Miller <davem at davemloft dot net>
- Date: Fri, 28 Sep 2012 20:05:40 +0000
- Subject: Fix sign of exact zero return from fma (bug 14638)
Bug 14638 reports fma producing a zero with the wrong sign for certain
cases where the result is an exact zero.
This patch fixes this by adding appropriate checks for the relevant
special cases. For some implementations, only the case where the
expression becomes 0 + 0 needs special handling; for other
implementations, the general case of adding two values with equal
magnitudes and opposite sign needs to be checked for.
Tests are added for the nondefault rounding modes; these are just a
starting point for testing fma in those rounding modes. I can see
cases in the implementations where there are likely to be bugs in
those rounding modes; it probably makes sense to fix those before
locating the latest version of Jakub's random test generation for fma
and getting it to test nondefault rounding modes so as to track down
other such bugs.
I also think there may be bugs in some of the implementations in the
sign of zero they return when it's an inexact zero (underflowing
multiplication, added to z == 0) but also consider that a separate
bug.
Tested x86_64 and x86. I'd welcome testing for an ldbl-128
architecture with rounding mode support, such as sparc, to verify that
I did find the right places to change for ldbl-128 (the ldbl-128
target most convenient for me to test on myself is mips64, which
doesn't have rounding modes or exceptions for long double).
2012-09-28 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.
diff --git a/math/libm-test.inc b/math/libm-test.inc
index e8398bd..007eea1 100644
--- 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);
@@ -4608,6 +4638,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)
{
START (fmax);
@@ -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 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);
diff --git a/sysdeps/ieee754/ldbl-128/s_fma.c b/sysdeps/ieee754/ldbl-128/s_fma.c
index 355b60e..b08ff29 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 963bbd7..df68ade 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 78c0b0d..001d806 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 ca1e090..c27b0bd 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;
--
Joseph S. Myers
joseph@codesourcery.com