Fix fma missing underflows and bad results for some subnormal results (bugs 14152, 14783).

3 files changed, 24 insertions, 42 deletions

diff --git a/sysdeps/ieee754/dbl-64/s_fma.c b/sysdeps/ieee754/dbl-64/s_fma.c
index 5e21461a4b..108fd66192 100644
--- a/sysdeps/ieee754/dbl-64/s_fma.c
+++ b/sysdeps/ieee754/dbl-64/s_fma.c
@@ -210,20 +210,14 @@ __fma (double x, double y, double z)
 	{
 	  /* v.ieee.mantissa1 & 2 is LSB bit of the result before rounding,
 	     v.ieee.mantissa1 & 1 is the round bit and j is our sticky
-	     bit.  In round-to-nearest 001 rounds down like 00,
-	     011 rounds up, even though 01 rounds down (thus we need
-	     to adjust), 101 rounds down like 10 and 111 rounds up
-	     like 11.  */
-	  if ((v.ieee.mantissa1 & 3) == 1)
-	    {
-	      v.d *= 0x1p-106;
-	      if (v.ieee.negative)
-		return v.d - 0x1p-1074 /* __DBL_DENORM_MIN__ */;
-	      else
-		return v.d + 0x1p-1074 /* __DBL_DENORM_MIN__ */;
-	    }
-	  else
-	    return v.d * 0x1p-106;
+	     bit.  */
+	  w.d = 0.0;
+	  w.ieee.mantissa1 = ((v.ieee.mantissa1 & 3) << 1) | j;
+	  w.ieee.negative = v.ieee.negative;
+	  v.ieee.mantissa1 &= ~3U;
+	  v.d *= 0x1p-106;
+	  w.d *= 0x1p-2;
+	  return v.d + w.d;
 	}
       v.ieee.mantissa1 |= j;
       return v.d * 0x1p-106;
diff --git a/sysdeps/ieee754/ldbl-128/s_fmal.c b/sysdeps/ieee754/ldbl-128/s_fmal.c
index 46b3d81ce5..7fb9856ef1 100644
--- a/sysdeps/ieee754/ldbl-128/s_fmal.c
+++ b/sysdeps/ieee754/ldbl-128/s_fmal.c
@@ -211,20 +211,14 @@ __fmal (long double x, long double y, long double z)
 	{
 	  /* v.ieee.mantissa3 & 2 is LSB bit of the result before rounding,
 	     v.ieee.mantissa3 & 1 is the round bit and j is our sticky
-	     bit.  In round-to-nearest 001 rounds down like 00,
-	     011 rounds up, even though 01 rounds down (thus we need
-	     to adjust), 101 rounds down like 10 and 111 rounds up
-	     like 11.  */
-	  if ((v.ieee.mantissa3 & 3) == 1)
-	    {
-	      v.d *= 0x1p-226L;
-	      if (v.ieee.negative)
-		return v.d - 0x1p-16494L /* __LDBL_DENORM_MIN__ */;
-	      else
-		return v.d + 0x1p-16494L /* __LDBL_DENORM_MIN__ */;
-	    }
-	  else
-	    return v.d * 0x1p-226L;
+	     bit.  */
+	  w.d = 0.0L;
+	  w.ieee.mantissa3 = ((v.ieee.mantissa3 & 3) << 1) | j;
+	  w.ieee.negative = v.ieee.negative;
+	  v.ieee.mantissa3 &= ~3U;
+	  v.d *= 0x1p-226L;
+	  w.d *= 0x1p-2L;
+	  return v.d + w.d;
 	}
       v.ieee.mantissa3 |= j;
       return v.d * 0x1p-226L;
diff --git a/sysdeps/ieee754/ldbl-96/s_fmal.c b/sysdeps/ieee754/ldbl-96/s_fmal.c
index d125124286..73104914b3 100644
--- a/sysdeps/ieee754/ldbl-96/s_fmal.c
+++ b/sysdeps/ieee754/ldbl-96/s_fmal.c
@@ -211,20 +211,14 @@ __fmal (long double x, long double y, long double z)
 	{
 	  /* v.ieee.mantissa1 & 2 is LSB bit of the result before rounding,
 	     v.ieee.mantissa1 & 1 is the round bit and j is our sticky
-	     bit.  In round-to-nearest 001 rounds down like 00,
-	     011 rounds up, even though 01 rounds down (thus we need
-	     to adjust), 101 rounds down like 10 and 111 rounds up
-	     like 11.  */
-	  if ((v.ieee.mantissa1 & 3) == 1)
-	    {
-	      v.d *= 0x1p-128L;
-	      if (v.ieee.negative)
-		return v.d - 0x1p-16445L /* __LDBL_DENORM_MIN__ */;
-	      else
-		return v.d + 0x1p-16445L /* __LDBL_DENORM_MIN__ */;
-	    }
-	  else
-	    return v.d * 0x1p-128L;
+	     bit.  */
+	  w.d = 0.0L;
+	  w.ieee.mantissa1 = ((v.ieee.mantissa1 & 3) << 1) | j;
+	  w.ieee.negative = v.ieee.negative;
+	  v.ieee.mantissa1 &= ~3U;
+	  v.d *= 0x1p-128L;
+	  w.d *= 0x1p-2L;
+	  return v.d + w.d;
 	}
       v.ieee.mantissa1 |= j;
       return v.d * 0x1p-128L;