Fix sign of inexact zero return from fma (bug 14645).

6 files changed, 74 insertions, 1 deletions

diff --git a/ChangeLog b/ChangeLog
index 9e8fca0b0f..4538135505 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,17 @@
+2012-10-01  Joseph Myers  <joseph@codesourcery.com>
+
+	[BZ #14645]
+	* sysdeps/ieee754/dbl-64/s_fma.c (__fma): Compute result as x * y
+	if x * y might underflow to zero and z is zero.
+	* sysdeps/ieee754/ldbl-128/s_fmal.c (__fmal): Likewise.
+	* sysdeps/ieee754/ldbl-96/s_fmal.c (__fmal): Likewise.
+	* math/libm-test.inc (min_subnorm_value): New variable.
+	(fma_test): Add more tests.
+	(fma_test_towardzero): Likewise.
+	(fma_test_downward): Likewise
+	(fma_test_upward): Likewise.
+	(initialize): Set min_subnorm_value.
+
 2012-09-29  Joseph Myers  <joseph@codesourcery.com>
 
 	[BZ #14638]
diff --git a/NEWS b/NEWS
index e816c23069..03a853a0de 100644
--- 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, 14638.
+  14583, 14587, 14621, 14638, 14645.
 
 * 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 007eea1f30..bed8fc6ecc 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -214,6 +214,7 @@ static int ignore_max_ulp;	/* Should we ignore max_ulp?  */
 
 static FLOAT minus_zero, plus_zero;
 static FLOAT plus_infty, minus_infty, nan_value, max_value, min_value;
+static FLOAT min_subnorm_value;
 
 static FLOAT max_error, real_max_error, imag_max_error;
 
@@ -4576,6 +4577,15 @@ fma_test (void)
   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, min_value, min_value, plus_zero, plus_zero, UNDERFLOW_EXCEPTION);
+  TEST_fff_f (fma, min_value, min_value, minus_zero, plus_zero, UNDERFLOW_EXCEPTION);
+  TEST_fff_f (fma, min_value, -min_value, plus_zero, minus_zero, UNDERFLOW_EXCEPTION);
+  TEST_fff_f (fma, min_value, -min_value, minus_zero, minus_zero, UNDERFLOW_EXCEPTION);
+  TEST_fff_f (fma, -min_value, min_value, plus_zero, minus_zero, UNDERFLOW_EXCEPTION);
+  TEST_fff_f (fma, -min_value, min_value, minus_zero, minus_zero, UNDERFLOW_EXCEPTION);
+  TEST_fff_f (fma, -min_value, -min_value, plus_zero, plus_zero, UNDERFLOW_EXCEPTION);
+  TEST_fff_f (fma, -min_value, -min_value, minus_zero, plus_zero, UNDERFLOW_EXCEPTION);
+
 #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);
@@ -4676,6 +4686,15 @@ fma_test_towardzero (void)
       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, min_value, min_value, plus_zero, plus_zero, UNDERFLOW_EXCEPTION);
+      TEST_fff_f (fma, min_value, min_value, minus_zero, plus_zero, UNDERFLOW_EXCEPTION);
+      TEST_fff_f (fma, min_value, -min_value, plus_zero, minus_zero, UNDERFLOW_EXCEPTION);
+      TEST_fff_f (fma, min_value, -min_value, minus_zero, minus_zero, UNDERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -min_value, min_value, plus_zero, minus_zero, UNDERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -min_value, min_value, minus_zero, minus_zero, UNDERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -min_value, -min_value, plus_zero, plus_zero, UNDERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -min_value, -min_value, minus_zero, plus_zero, UNDERFLOW_EXCEPTION);
     }
 
   fesetround (save_round_mode);
@@ -4723,6 +4742,15 @@ fma_test_downward (void)
       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, min_value, min_value, plus_zero, plus_zero, UNDERFLOW_EXCEPTION);
+      TEST_fff_f (fma, min_value, min_value, minus_zero, plus_zero, UNDERFLOW_EXCEPTION);
+      TEST_fff_f (fma, min_value, -min_value, plus_zero, -min_subnorm_value, UNDERFLOW_EXCEPTION);
+      TEST_fff_f (fma, min_value, -min_value, minus_zero, -min_subnorm_value, UNDERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -min_value, min_value, plus_zero, -min_subnorm_value, UNDERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -min_value, min_value, minus_zero, -min_subnorm_value, UNDERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -min_value, -min_value, plus_zero, plus_zero, UNDERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -min_value, -min_value, minus_zero, plus_zero, UNDERFLOW_EXCEPTION);
     }
 
   fesetround (save_round_mode);
@@ -4770,6 +4798,15 @@ fma_test_upward (void)
       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, min_value, min_value, plus_zero, min_subnorm_value, UNDERFLOW_EXCEPTION);
+      TEST_fff_f (fma, min_value, min_value, minus_zero, min_subnorm_value, UNDERFLOW_EXCEPTION);
+      TEST_fff_f (fma, min_value, -min_value, plus_zero, minus_zero, UNDERFLOW_EXCEPTION);
+      TEST_fff_f (fma, min_value, -min_value, minus_zero, minus_zero, UNDERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -min_value, min_value, plus_zero, minus_zero, UNDERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -min_value, min_value, minus_zero, minus_zero, UNDERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -min_value, -min_value, plus_zero, min_subnorm_value, UNDERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -min_value, -min_value, minus_zero, min_subnorm_value, UNDERFLOW_EXCEPTION);
     }
 
   fesetround (save_round_mode);
@@ -9431,6 +9468,12 @@ initialize (void)
 		      LDBL_MAX, DBL_MAX, FLT_MAX);
   min_value = CHOOSE (LDBL_MIN, DBL_MIN, FLT_MIN,
 		      LDBL_MIN, DBL_MIN, FLT_MIN);
+  min_subnorm_value = CHOOSE (__LDBL_DENORM_MIN__,
+			      __DBL_DENORM_MIN__,
+			      __FLT_DENORM_MIN__,
+			      __LDBL_DENORM_MIN__,
+			      __DBL_DENORM_MIN__,
+			      __FLT_DENORM_MIN__);
 
   (void) &plus_zero;
   (void) &nan_value;
@@ -9439,6 +9482,7 @@ initialize (void)
   (void) &minus_infty;
   (void) &max_value;
   (void) &min_value;
+  (void) &min_subnorm_value;
 
   /* Clear all exceptions.  From now on we must not get random exceptions.  */
   feclearexcept (FE_ALL_EXCEPT);
diff --git a/sysdeps/ieee754/dbl-64/s_fma.c b/sysdeps/ieee754/dbl-64/s_fma.c
index c9809fb180..5e21461a4b 100644
--- a/sysdeps/ieee754/dbl-64/s_fma.c
+++ b/sysdeps/ieee754/dbl-64/s_fma.c
@@ -49,6 +49,11 @@ __fma (double x, double y, double z)
 	  && u.ieee.exponent != 0x7ff
 	  && v.ieee.exponent != 0x7ff)
 	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 DBL_DENORM_MIN,
 	 compute as x * y + z.  */
diff --git a/sysdeps/ieee754/ldbl-128/s_fmal.c b/sysdeps/ieee754/ldbl-128/s_fmal.c
index df68ade435..46b3d81ce5 100644
--- a/sysdeps/ieee754/ldbl-128/s_fmal.c
+++ b/sysdeps/ieee754/ldbl-128/s_fmal.c
@@ -50,6 +50,11 @@ __fmal (long double x, long double y, long double 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 LDBL_DENORM_MIN,
 	 compute as x * y + z.  */
diff --git a/sysdeps/ieee754/ldbl-96/s_fmal.c b/sysdeps/ieee754/ldbl-96/s_fmal.c
index c27b0bd852..d125124286 100644
--- a/sysdeps/ieee754/ldbl-96/s_fmal.c
+++ b/sysdeps/ieee754/ldbl-96/s_fmal.c
@@ -50,6 +50,11 @@ __fmal (long double x, long double y, long double 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 LDBL_DENORM_MIN,
 	 compute as x * y + z.  */