[PATCH 3/4] math: signal underflow for fmaf results that are tiny before rounding

fmaf computes the exact result as a double and lets the return convert
it to float, so that conversion is the operation's single rounding and
is what has to signal underflow.  Where the exact result is tiny but
rounds up to the smallest normal float, architectures determining
tininess after rounding do not signal it, because the value they examine
is no longer tiny.

Route the returns through a helper that applies CHECK_NARROW_TINY, which
already decides this by rounding with an unbounded exponent range.  The
raise is redundant on architectures that signal such results themselves.
---
 sysdeps/ieee754/dbl-64/s_fmaf.c | 23 +++++++++++++++++++----
 1 file changed, 19 insertions(+), 4 deletions(-)

diff --git ./sysdeps/ieee754/dbl-64/s_fmaf.c ./sysdeps/ieee754/dbl-64/s_fmaf.c
index 639c273486..dc78e9453a 100644
--- ./sysdeps/ieee754/dbl-64/s_fmaf.c
+++ ./sysdeps/ieee754/dbl-64/s_fmaf.c
@@ -20,9 +20,24 @@
 #include <math.h>
 #include <fenv.h>
 #include <libm-alias-float.h>
+#include <math-narrow.h>
 #include <math-use-builtins.h>
 #include "math_config.h"
 
+#if !USE_FMAF_BUILTIN
+/* Narrow the double result to float.  The exact result of the fma is
+   representable as a double, so the conversion is the operation's single
+   rounding and it is what must signal underflow; help it where the
+   architecture does not.  */
+static inline float
+narrow_fmaf_result (double d)
+{
+  float ret = (float) d;
+  CHECK_NARROW_TINY (ret, d, __FLT_MIN__);
+  return ret;
+}
+#endif
+
 float
 __fmaf (float x, float y, float z)
 {
@@ -37,17 +52,17 @@ __fmaf (float x, float y, float z)
   /* If not exact or at round to even boundary, the result is correct in
      all rounding modes.  */
   if (__glibc_likely ((u & 0xfffffff) != 0))
-    return result;
+    return narrow_fmaf_result (result);
 
   /* Also check if the double result appears exact when it might not be and
      thus it will not set the underflow flag if denormal.  */
   if ((u & 0x10000000) == 0
       && ((u >> MANTISSA_WIDTH) & 0x7ff) > EXPONENT_BIAS - 126)
-    return result;
+    return narrow_fmaf_result (result);
 
   /* Return if result is exact in all rounding modes.  */
   if (result - xy == z && result - z == xy)
-    return result;
+    return narrow_fmaf_result (result);
 
   /* This is where 'double-rouding' might return a wrong value, and thus
      needs adjusting the low-order bits in the direction of the error.  */
@@ -62,7 +77,7 @@ __fmaf (float x, float y, float z)
     u++;
   else
     u--;
-  return asdouble (u);
+  return narrow_fmaf_result (asdouble (u));
 #endif /* ! USE_FMAF_BUILTIN  */
 }
 #ifndef __fmaf
-- 
2.54.0