Fix sign of inexact zero return from fma (bug 14645)
Joseph S. Myers
joseph@codesourcery.com
Sat Sep 29 21:51:00 GMT 2012
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
More information about the Libc-alpha
mailing list