Fix fma overflow results outside round-to-nearest mode (bug 14797).

diff --git a/ChangeLog b/ChangeLog
index a1a7a04927be86914c1bc66c2bc013280918ad9d..9b5b8f79bafdeb13948266ffa1eb41e8f1f9e372 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,15 @@
+2012-11-04  Joseph Myers  <joseph@codesourcery.com>
+
+       [BZ #14797]
+       * sysdeps/ieee754/dbl-64/s_fma.c (__fma): Compute cases that
+       definitely overflow as x * y not x * y + z.
+       * sysdeps/ieee754/ldbl-128/s_fmal.c (__fmal): Likewise.
+       * sysdeps/ieee754/ldbl-96/s_fmal.c (__fmal): Likewise.
+       * math/libm-test.inc (fma_test): Add more tests.
+       (fma_test_towardzero): Likewise.
+       (fma_test_downward): Likewise.
+       (fma_test_upward): Likewise.
+
 2012-11-04  Thomas Schwinge  <thomas@codesourcery.com>
 
        [BZ #157]

diff --git a/NEWS b/NEWS
index e823c2bcf7c4a51ae59b4b15316ad0a1a885b25b..a6d757f153bc6c9e5646ff7da163948f952d7257 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -18,7 +18,7 @@ Version 2.17
   14530, 14532, 14538, 14543, 14544, 14545, 14557, 14562, 14568, 14576,
   14579, 14583, 14587, 14595, 14602, 14610, 14621, 14638, 14645, 14648,
   14652, 14660, 14661, 14669, 14683, 14694, 14716, 14743, 14767, 14783,
-  14784, 14785, 14796.
+  14784, 14785, 14796, 14797.
 
 * 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 48241e02f899cfba8c14a106331423536774d55f..55892c3459bb908fbb24f5fafafb6a6713c8a087 100644 (file)
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -4617,6 +4617,15 @@ fma_test (void)
   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, max_value, max_value, min_value, plus_infty, OVERFLOW_EXCEPTION);
+  TEST_fff_f (fma, max_value, max_value, -min_value, plus_infty, OVERFLOW_EXCEPTION);
+  TEST_fff_f (fma, max_value, -max_value, min_value, minus_infty, OVERFLOW_EXCEPTION);
+  TEST_fff_f (fma, max_value, -max_value, -min_value, minus_infty, OVERFLOW_EXCEPTION);
+  TEST_fff_f (fma, -max_value, max_value, min_value, minus_infty, OVERFLOW_EXCEPTION);
+  TEST_fff_f (fma, -max_value, max_value, -min_value, minus_infty, OVERFLOW_EXCEPTION);
+  TEST_fff_f (fma, -max_value, -max_value, min_value, plus_infty, OVERFLOW_EXCEPTION);
+  TEST_fff_f (fma, -max_value, -max_value, -min_value, plus_infty, OVERFLOW_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);
@@ -4837,6 +4846,15 @@ fma_test_towardzero (void)
       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, max_value, max_value, min_value, max_value, OVERFLOW_EXCEPTION);
+      TEST_fff_f (fma, max_value, max_value, -min_value, max_value, OVERFLOW_EXCEPTION);
+      TEST_fff_f (fma, max_value, -max_value, min_value, -max_value, OVERFLOW_EXCEPTION);
+      TEST_fff_f (fma, max_value, -max_value, -min_value, -max_value, OVERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -max_value, max_value, min_value, -max_value, OVERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -max_value, max_value, -min_value, -max_value, OVERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -max_value, -max_value, min_value, max_value, OVERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -max_value, -max_value, -min_value, max_value, OVERFLOW_EXCEPTION);
+
 #if defined (TEST_FLOAT) && FLT_MANT_DIG == 24
       TEST_fff_f (fma, 0x1.4p-126, 0x1.000004p-1, 0x1p-128, 0x1.c00004p-127, UNDERFLOW_EXCEPTION);
       TEST_fff_f (fma, -0x1.4p-126, 0x1.000004p-1, -0x1p-128, -0x1.c00004p-127, UNDERFLOW_EXCEPTION);
@@ -5014,6 +5032,15 @@ fma_test_downward (void)
       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, max_value, max_value, min_value, max_value, OVERFLOW_EXCEPTION);
+      TEST_fff_f (fma, max_value, max_value, -min_value, max_value, OVERFLOW_EXCEPTION);
+      TEST_fff_f (fma, max_value, -max_value, min_value, minus_infty, OVERFLOW_EXCEPTION);
+      TEST_fff_f (fma, max_value, -max_value, -min_value, minus_infty, OVERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -max_value, max_value, min_value, minus_infty, OVERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -max_value, max_value, -min_value, minus_infty, OVERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -max_value, -max_value, min_value, max_value, OVERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -max_value, -max_value, -min_value, max_value, OVERFLOW_EXCEPTION);
+
 #if defined (TEST_FLOAT) && FLT_MANT_DIG == 24
       TEST_fff_f (fma, 0x1.4p-126, 0x1.000004p-1, 0x1p-128, 0x1.c00004p-127, UNDERFLOW_EXCEPTION);
       TEST_fff_f (fma, -0x1.4p-126, 0x1.000004p-1, -0x1p-128, -0x1.c00008p-127, UNDERFLOW_EXCEPTION);
@@ -5191,6 +5218,15 @@ fma_test_upward (void)
       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, max_value, max_value, min_value, plus_infty, OVERFLOW_EXCEPTION);
+      TEST_fff_f (fma, max_value, max_value, -min_value, plus_infty, OVERFLOW_EXCEPTION);
+      TEST_fff_f (fma, max_value, -max_value, min_value, -max_value, OVERFLOW_EXCEPTION);
+      TEST_fff_f (fma, max_value, -max_value, -min_value, -max_value, OVERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -max_value, max_value, min_value, -max_value, OVERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -max_value, max_value, -min_value, -max_value, OVERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -max_value, -max_value, min_value, plus_infty, OVERFLOW_EXCEPTION);
+      TEST_fff_f (fma, -max_value, -max_value, -min_value, plus_infty, OVERFLOW_EXCEPTION);
+
 #if defined (TEST_FLOAT) && FLT_MANT_DIG == 24
       TEST_fff_f (fma, 0x1.4p-126, 0x1.000004p-1, 0x1p-128, 0x1.c00008p-127, UNDERFLOW_EXCEPTION);
       TEST_fff_f (fma, -0x1.4p-126, 0x1.000004p-1, -0x1p-128, -0x1.c00004p-127, UNDERFLOW_EXCEPTION);

diff --git a/sysdeps/ieee754/dbl-64/s_fma.c b/sysdeps/ieee754/dbl-64/s_fma.c
index 07fe715617f69fc76df0f1857d5c29ad5b4f3b3b..cd2883070950fb970a5e9f76aaf19e5eb7594991 100644 (file)
--- a/sysdeps/ieee754/dbl-64/s_fma.c
+++ b/sysdeps/ieee754/dbl-64/s_fma.c
@@ -55,16 +55,17 @@ __fma (double x, double y, double z)
         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 zero, compute as x * y + z.  */
+      /* If x or y or z is Inf/NaN, or if x * y is zero, compute as
+        x * y + z.  */
       if (u.ieee.exponent == 0x7ff
          || v.ieee.exponent == 0x7ff
          || w.ieee.exponent == 0x7ff
-         || u.ieee.exponent + v.ieee.exponent
-            > 0x7ff + IEEE754_DOUBLE_BIAS
          || x == 0
          || y == 0)
        return x * y + z;
+      /* If fma will certainly overflow, compute as x * y.  */
+      if (u.ieee.exponent + v.ieee.exponent > 0x7ff + IEEE754_DOUBLE_BIAS)
+       return x * y;
       /* If x * y is less than 1/4 of DBL_DENORM_MIN, neither the
         result nor whether there is underflow depends on its exact
         value, only on its sign.  */

diff --git a/sysdeps/ieee754/ldbl-128/s_fmal.c b/sysdeps/ieee754/ldbl-128/s_fmal.c
index d576403b2582df2e5f022342817325dc5d6f918a..6fa663a6c09e94c6c84a35d3ba81bcc8c3c4c5ab 100644 (file)
--- a/sysdeps/ieee754/ldbl-128/s_fmal.c
+++ b/sysdeps/ieee754/ldbl-128/s_fmal.c
@@ -56,16 +56,18 @@ __fmal (long double x, long double y, long double z)
         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 zero, compute as x * y + z.  */
+      /* If x or y or z is Inf/NaN, or if x * y is zero, compute as
+        x * y + z.  */
       if (u.ieee.exponent == 0x7fff
          || v.ieee.exponent == 0x7fff
          || w.ieee.exponent == 0x7fff
-         || u.ieee.exponent + v.ieee.exponent
-            > 0x7fff + IEEE854_LONG_DOUBLE_BIAS
          || x == 0
          || y == 0)
        return x * y + z;
+      /* If fma will certainly overflow, compute as x * y.  */
+      if (u.ieee.exponent + v.ieee.exponent
+         > 0x7fff + IEEE854_LONG_DOUBLE_BIAS)
+       return x * y;
       /* If x * y is less than 1/4 of LDBL_DENORM_MIN, neither the
         result nor whether there is underflow depends on its exact
         value, only on its sign.  */

diff --git a/sysdeps/ieee754/ldbl-96/s_fmal.c b/sysdeps/ieee754/ldbl-96/s_fmal.c
index 32e71a18ba966f6f21b573f5598471aa0c5ce048..53098b6d4efb562ec08cb49d7643d5baa16c9f91 100644 (file)
--- a/sysdeps/ieee754/ldbl-96/s_fmal.c
+++ b/sysdeps/ieee754/ldbl-96/s_fmal.c
@@ -56,16 +56,18 @@ __fmal (long double x, long double y, long double z)
         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 zero, compute as x * y + z.  */
+      /* If x or y or z is Inf/NaN, or if x * y is zero, compute as
+        x * y + z.  */
       if (u.ieee.exponent == 0x7fff
          || v.ieee.exponent == 0x7fff
          || w.ieee.exponent == 0x7fff
-         || u.ieee.exponent + v.ieee.exponent
-            > 0x7fff + IEEE854_LONG_DOUBLE_BIAS
          || x == 0
          || y == 0)
        return x * y + z;
+      /* If fma will certainly overflow, compute as x * y.  */
+      if (u.ieee.exponent + v.ieee.exponent
+         > 0x7fff + IEEE854_LONG_DOUBLE_BIAS)
+       return x * y;
       /* If x * y is less than 1/4 of LDBL_DENORM_MIN, neither the
         result nor whether there is underflow depends on its exact
         value, only on its sign.  */