This is the mail archive of the libc-alpha@sourceware.org mailing list for the glibc project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

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


Bug 14645 is incorrect sign of inexact zero results from fma
functions, where the multiplication underflows to zero so the code
assumes that x * y + z is an appropriate computation, but because z is
also zero the addition ends up with zero of a sign determined by the
rounding-mode-dependent rules for adding 0 + 0, when the result should
always have the sign of the infinite-precision value of x * y in this
case.

This patch adds an appropriate check for this case and uses x * y
instead of x * y + z as the underlying computation.  (I think there
are still related bugs involving spurious exceptions, and wrong
results in round-to-zero mode, from the use of x * y + z when x * y
underflows and z is finite but nonzero; again I consider those
separate bugs.)

Tested x86_64 (both multi-arch enabled, and multi-arch disabled) and
x86.  Again, I'd welcome testing for an ldbl-128 architecture with
rounding mode support, such as sparc, to confirm that the changes are
sufficient for ldbl-128 (for these particular tests).

2012-09-29  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.

diff --git a/math/libm-test.inc b/math/libm-test.inc
index 007eea1..bed8fc6 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 c9809fb..5e21461 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 df68ade..46b3d81 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 c27b0bd..d125124 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.  */

-- 
Joseph S. Myers
joseph@codesourcery.com


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]