Fix f64xdivf128, f64xmulf128 spurious underflows (bug 28358)

As described in bug 28358, the round-to-odd computations used in the
libm functions that round their results to a narrower format can yield
spurious underflow exceptions in the following circumstances: the
narrowing only narrows the precision of the type and not the exponent
range (i.e., it's narrowing _Float128 to _Float64x on x86_64, x86 or
ia64), the architecture does after-rounding tininess detection (which
applies to all those architectures), the result is inexact, tiny
before rounding but not tiny after rounding (with the chosen rounding
mode) for _Float64x (which is possible for narrowing mul, div and fma,
not for narrowing add, sub or sqrt), so the underflow exception
resulting from the toward-zero computation in _Float128 is spurious
for _Float64x.

Fixed by making ROUND_TO_ODD call feclearexcept (FE_UNDERFLOW) in the
problem cases (as indicated by an extra argument to the macro); there
is never any need to preserve underflow exceptions from this part of
the computation, because the conversion of the round-to-odd value to
the narrower type will underflow in exactly the cases in which the
function should raise that exception, but it may be more efficient to
avoid the extra manipulation of the floating-point environment when
not needed.

Tested for x86_64 and x86, and with build-many-glibcs.py.

12 files changed, 16 insertions, 12 deletions

diff --git a/sysdeps/ieee754/dbl-64/s_fdiv.c b/sysdeps/ieee754/dbl-64/s_fdiv.c
index e433aee309..98640a235c 100644
--- a/sysdeps/ieee754/dbl-64/s_fdiv.c
+++ b/sysdeps/ieee754/dbl-64/s_fdiv.c
@@ -29,6 +29,7 @@
 float
 __fdiv (double x, double y)
 {
-  NARROW_DIV_ROUND_TO_ODD (x, y, float, union ieee754_double, , mantissa1);
+  NARROW_DIV_ROUND_TO_ODD (x, y, float, union ieee754_double, , mantissa1,
+			   false);
 }
 libm_alias_float_double (div)
diff --git a/sysdeps/ieee754/dbl-64/s_fmul.c b/sysdeps/ieee754/dbl-64/s_fmul.c
index 892776f2d3..d97bbf60ed 100644
--- a/sysdeps/ieee754/dbl-64/s_fmul.c
+++ b/sysdeps/ieee754/dbl-64/s_fmul.c
@@ -29,6 +29,7 @@
 float
 __fmul (double x, double y)
 {
-  NARROW_MUL_ROUND_TO_ODD (x, y, float, union ieee754_double, , mantissa1);
+  NARROW_MUL_ROUND_TO_ODD (x, y, float, union ieee754_double, , mantissa1,
+			   false);
 }
 libm_alias_float_double (mul)
diff --git a/sysdeps/ieee754/ldbl-128/s_ddivl.c b/sysdeps/ieee754/ldbl-128/s_ddivl.c
index 80c58f07e9..88e242bccb 100644
--- a/sysdeps/ieee754/ldbl-128/s_ddivl.c
+++ b/sysdeps/ieee754/ldbl-128/s_ddivl.c
@@ -32,6 +32,6 @@ double
 __ddivl (_Float128 x, _Float128 y)
 {
   NARROW_DIV_ROUND_TO_ODD (x, y, double, union ieee854_long_double, l,
-			   mantissa3);
+			   mantissa3, false);
 }
 libm_alias_double_ldouble (div)
diff --git a/sysdeps/ieee754/ldbl-128/s_dmull.c b/sysdeps/ieee754/ldbl-128/s_dmull.c
index 9b71bbab1e..186ba172bb 100644
--- a/sysdeps/ieee754/ldbl-128/s_dmull.c
+++ b/sysdeps/ieee754/ldbl-128/s_dmull.c
@@ -32,6 +32,6 @@ double
 __dmull (_Float128 x, _Float128 y)
 {
   NARROW_MUL_ROUND_TO_ODD (x, y, double, union ieee854_long_double, l,
-			   mantissa3);
+			   mantissa3, false);
 }
 libm_alias_double_ldouble (mul)
diff --git a/sysdeps/ieee754/ldbl-128/s_f64xdivf128.c b/sysdeps/ieee754/ldbl-128/s_f64xdivf128.c
index a6f114c01c..03f2052adf 100644
--- a/sysdeps/ieee754/ldbl-128/s_f64xdivf128.c
+++ b/sysdeps/ieee754/ldbl-128/s_f64xdivf128.c
@@ -18,6 +18,7 @@
 
 #include <math.h>
 #include <math-narrow.h>
+#include <tininess.h>
 
 /* math_ldbl.h defines _Float128 to long double for this directory,
    but when they are different, this function must be defined with
@@ -30,7 +31,7 @@ __f64xdivf128 (_Float128 x, _Float128 y)
 {
 #if __HAVE_FLOAT64X_LONG_DOUBLE && __HAVE_DISTINCT_FLOAT128
   NARROW_DIV_ROUND_TO_ODD (x, y, _Float64x, union ieee854_long_double, l,
-			   mantissa3);
+			   mantissa3, TININESS_AFTER_ROUNDING);
 #else
   NARROW_DIV_TRIVIAL (x, y, _Float64x);
 #endif
diff --git a/sysdeps/ieee754/ldbl-128/s_f64xmulf128.c b/sysdeps/ieee754/ldbl-128/s_f64xmulf128.c
index ebb34ad4e0..031f8ef480 100644
--- a/sysdeps/ieee754/ldbl-128/s_f64xmulf128.c
+++ b/sysdeps/ieee754/ldbl-128/s_f64xmulf128.c
@@ -18,6 +18,7 @@
 
 #include <math.h>
 #include <math-narrow.h>
+#include <tininess.h>
 
 /* math_ldbl.h defines _Float128 to long double for this directory,
    but when they are different, this function must be defined with
@@ -30,7 +31,7 @@ __f64xmulf128 (_Float128 x, _Float128 y)
 {
 #if __HAVE_FLOAT64X_LONG_DOUBLE && __HAVE_DISTINCT_FLOAT128
   NARROW_MUL_ROUND_TO_ODD (x, y, _Float64x, union ieee854_long_double, l,
-			   mantissa3);
+			   mantissa3, TININESS_AFTER_ROUNDING);
 #else
   NARROW_MUL_TRIVIAL (x, y, _Float64x);
 #endif
diff --git a/sysdeps/ieee754/ldbl-128/s_fdivl.c b/sysdeps/ieee754/ldbl-128/s_fdivl.c
index dbad331909..7ceb4db4f8 100644
--- a/sysdeps/ieee754/ldbl-128/s_fdivl.c
+++ b/sysdeps/ieee754/ldbl-128/s_fdivl.c
@@ -28,6 +28,6 @@ float
 __fdivl (_Float128 x, _Float128 y)
 {
   NARROW_DIV_ROUND_TO_ODD (x, y, float, union ieee854_long_double, l,
-			   mantissa3);
+			   mantissa3, false);
 }
 libm_alias_float_ldouble (div)
diff --git a/sysdeps/ieee754/ldbl-128/s_fmull.c b/sysdeps/ieee754/ldbl-128/s_fmull.c
index 1b6cb25921..b2602fcd61 100644
--- a/sysdeps/ieee754/ldbl-128/s_fmull.c
+++ b/sysdeps/ieee754/ldbl-128/s_fmull.c
@@ -28,6 +28,6 @@ float
 __fmull (_Float128 x, _Float128 y)
 {
   NARROW_MUL_ROUND_TO_ODD (x, y, float, union ieee854_long_double, l,
-			   mantissa3);
+			   mantissa3, false);
 }
 libm_alias_float_ldouble (mul)
diff --git a/sysdeps/ieee754/ldbl-96/s_ddivl.c b/sysdeps/ieee754/ldbl-96/s_ddivl.c
index 5ec7a78d8a..32df44da24 100644
--- a/sysdeps/ieee754/ldbl-96/s_ddivl.c
+++ b/sysdeps/ieee754/ldbl-96/s_ddivl.c
@@ -28,6 +28,6 @@ double
 __ddivl (long double x, long double y)
 {
   NARROW_DIV_ROUND_TO_ODD (x, y, double, union ieee854_long_double, l,
-			   mantissa1);
+			   mantissa1, false);
 }
 libm_alias_double_ldouble (div)
diff --git a/sysdeps/ieee754/ldbl-96/s_dmull.c b/sysdeps/ieee754/ldbl-96/s_dmull.c
index b06bdfed52..b5bfb3f37b 100644
--- a/sysdeps/ieee754/ldbl-96/s_dmull.c
+++ b/sysdeps/ieee754/ldbl-96/s_dmull.c
@@ -28,6 +28,6 @@ double
 __dmull (long double x, long double y)
 {
   NARROW_MUL_ROUND_TO_ODD (x, y, double, union ieee854_long_double, l,
-			   mantissa1);
+			   mantissa1, false);
 }
 libm_alias_double_ldouble (mul)
diff --git a/sysdeps/ieee754/ldbl-96/s_fdivl.c b/sysdeps/ieee754/ldbl-96/s_fdivl.c
index 0b30cbb30a..9568b148bb 100644
--- a/sysdeps/ieee754/ldbl-96/s_fdivl.c
+++ b/sysdeps/ieee754/ldbl-96/s_fdivl.c
@@ -26,6 +26,6 @@ float
 __fdivl (long double x, long double y)
 {
   NARROW_DIV_ROUND_TO_ODD (x, y, float, union ieee854_long_double, l,
-			   mantissa1);
+			   mantissa1, false);
 }
 libm_alias_float_ldouble (div)
diff --git a/sysdeps/ieee754/ldbl-96/s_fmull.c b/sysdeps/ieee754/ldbl-96/s_fmull.c
index 6c9c9bc899..c3a5e1feb0 100644
--- a/sysdeps/ieee754/ldbl-96/s_fmull.c
+++ b/sysdeps/ieee754/ldbl-96/s_fmull.c
@@ -26,6 +26,6 @@ float
 __fmull (long double x, long double y)
 {
   NARROW_MUL_ROUND_TO_ODD (x, y, float, union ieee854_long_double, l,
-			   mantissa1);
+			   mantissa1, false);
 }
 libm_alias_float_ldouble (mul)