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

diff --git a/ChangeLog b/ChangeLog
index 9e8fca0b0f3d1396266542c9f71515c91a28bf75..453813550556ceb469b45fa4eee949b303577f85 100644 (file)
--- 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 e816c23069307493766d5eacaa345ba5692cf6be..03a853a0deb337fc9765daa0370ef24206a82da6 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, 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 007eea1f30658f10d440bda8ea3697e44ef3fa9c..bed8fc6ecc0542ccf8ed82ab24279d965a884759 100644 (file)
--- 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 c9809fb180685a941044f4907da4a275f69e3c6d..5e21461a4b6ba428434e8137b50e1c6fab4d5ecd 100644 (file)
--- 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 df68ade435f2707856c0b7e0cca8ffb8e6e4dd7d..46b3d81ce53b8bbc56e4381df79884c654c3e37a 100644 (file)
--- 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 c27b0bd852c92ee7de6212ddbfa4551724373c03..d1251242867b17490e15bc9747c8ce7101b65f00 100644 (file)
--- 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.  */