Fix fma spurious underflows (bug 18824) [committed]

[Joseph Myers <joseph at codesourcery dot com>]

Various fma implementations have logic that, when computing fma (x, y,
z) where z is large (so care needs taking to avoid internal overflow)
but x * y is small, scale x * y up instead of down to avoid internal
underflows resulting from scaling down.  (In these cases, x * y is
small enough that only its sign actually matters rather than the exact
value.)

The threshold for scaling up instead of down was correct for "if the
unscaled values were multiplied, the low part of the multiplication
could underflow", and the scaling was sufficient to ensure that the
low part of the multiplication did not underflow (given that cases of
very small x * y - less than half the least subnormal - were
previously dealt with).  However, the choice in the functions wasn't
between scaling up or no scaling, but between scaling up and scaling
down (scaling down actually being needed when x * y isn't so small
compared to z and so the exact value does matter).  Thus a larger
threshold is needed to ensure that scaling down doesn't produce values
the multiplication of whose low parts underflows.  This patch
increases the thresholds accordingly.

Tested for x86_64, x86 and mips64 (with the MIPS version of s_fmal.c
removed so that the ldbl-128 version gets tested instead of the
soft-fp one).  Committed.

(auto-libm-test-out diffs omitted below.)

2015-08-14  Joseph Myers  <joseph@codesourcery.com>

	[BZ #18824]
	* sysdeps/ieee754/dbl-64/s_fma.c (__fma): Increase threshold for
	scaling x * y up instead of down.
	* sysdeps/ieee754/ldbl-128/s_fmal.c (__fmal): Likewise.
	* sysdeps/ieee754/ldbl-96/s_fmal.c (__fmal): Likewise.
	* math/auto-libm-test-in: Add more tests of fma.
	* math/auto-libm-test-out: Regenerated.

diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index d3a4834..4480371 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -1794,6 +1794,24 @@ fma 0x1.0000000000000000000000000001p-16382 -0x1.0000000000000000000000000001p-6
 fma 0x1.0000000000000000000000000001p-16382 0x1.0000000000000000000000000001p-66 -0x1p16319
 fma 0x1.0000000000000000000000000001p-16382 -0x1.0000000000000000000000000001p-66 -0x1p16319
 
+# Bug 6801: errno setting may be missing.
+fma 0x1.fffffep-126 0x1.fffffep25 0x1.fffffep127 missing-errno
+fma 0x1.fffffep-126 -0x1.fffffep25 0x1.fffffep127
+fma 0x1.fffffep-126 0x1.fffffep25 -0x1.fffffep127
+fma 0x1.fffffep-126 -0x1.fffffep25 -0x1.fffffep127 missing-errno
+fma 0x1.fffffffffffffp-1022 0x1.fffffffffffffp54 0x1.fffffffffffffp1023 missing-errno
+fma 0x1.fffffffffffffp-1022 -0x1.fffffffffffffp54 0x1.fffffffffffffp1023
+fma 0x1.fffffffffffffp-1022 0x1.fffffffffffffp54 -0x1.fffffffffffffp1023
+fma 0x1.fffffffffffffp-1022 -0x1.fffffffffffffp54 -0x1.fffffffffffffp1023 missing-errno
+fma 0x1.fffffffffffffffep-16382 0x1.fffffffffffffffep65 0x1.fffffffffffffffep16383 missing-errno
+fma 0x1.fffffffffffffffep-16382 -0x1.fffffffffffffffep65 0x1.fffffffffffffffep16383
+fma 0x1.fffffffffffffffep-16382 0x1.fffffffffffffffep65 -0x1.fffffffffffffffep16383
+fma 0x1.fffffffffffffffep-16382 -0x1.fffffffffffffffep65 -0x1.fffffffffffffffep16383 missing-errno
+fma 0x1.ffffffffffffffffffffffffffffp-16382 0x1.ffffffffffffffffffffffffffffp114 0x1.ffffffffffffffffffffffffffffp16383 missing-errno
+fma 0x1.ffffffffffffffffffffffffffffp-16382 -0x1.ffffffffffffffffffffffffffffp114 0x1.ffffffffffffffffffffffffffffp16383
+fma 0x1.ffffffffffffffffffffffffffffp-16382 0x1.ffffffffffffffffffffffffffffp114 -0x1.ffffffffffffffffffffffffffffp16383
+fma 0x1.ffffffffffffffffffffffffffffp-16382 -0x1.ffffffffffffffffffffffffffffp114 -0x1.ffffffffffffffffffffffffffffp16383 missing-errno
+
 hypot 0 0
 hypot 0 -0
 hypot -0 0
diff --git a/sysdeps/ieee754/dbl-64/s_fma.c b/sysdeps/ieee754/dbl-64/s_fma.c
index 716b412..278b690 100644
--- a/sysdeps/ieee754/dbl-64/s_fma.c
+++ b/sysdeps/ieee754/dbl-64/s_fma.c
@@ -117,7 +117,7 @@ __fma (double x, double y, double z)
 	     very small, adjust them up to avoid spurious underflows,
 	     rather than down.  */
 	  if (u.ieee.exponent + v.ieee.exponent
-	      <= IEEE754_DOUBLE_BIAS + DBL_MANT_DIG)
+	      <= IEEE754_DOUBLE_BIAS + 2 * DBL_MANT_DIG)
 	    {
 	      if (u.ieee.exponent > v.ieee.exponent)
 		u.ieee.exponent += 2 * DBL_MANT_DIG + 2;
diff --git a/sysdeps/ieee754/ldbl-128/s_fmal.c b/sysdeps/ieee754/ldbl-128/s_fmal.c
index b13178f..5abc910 100644
--- a/sysdeps/ieee754/ldbl-128/s_fmal.c
+++ b/sysdeps/ieee754/ldbl-128/s_fmal.c
@@ -121,7 +121,7 @@ __fmal (long double x, long double y, long double z)
 	     very small, adjust them up to avoid spurious underflows,
 	     rather than down.  */
 	  if (u.ieee.exponent + v.ieee.exponent
-	      <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG)
+	      <= IEEE854_LONG_DOUBLE_BIAS + 2 * LDBL_MANT_DIG)
 	    {
 	      if (u.ieee.exponent > v.ieee.exponent)
 		u.ieee.exponent += 2 * LDBL_MANT_DIG + 2;
diff --git a/sysdeps/ieee754/ldbl-96/s_fmal.c b/sysdeps/ieee754/ldbl-96/s_fmal.c
index eec5a02..1232c9e 100644
--- a/sysdeps/ieee754/ldbl-96/s_fmal.c
+++ b/sysdeps/ieee754/ldbl-96/s_fmal.c
@@ -119,7 +119,7 @@ __fmal (long double x, long double y, long double z)
 	     very small, adjust them up to avoid spurious underflows,
 	     rather than down.  */
 	  if (u.ieee.exponent + v.ieee.exponent
-	      <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG)
+	      <= IEEE854_LONG_DOUBLE_BIAS + 2 * LDBL_MANT_DIG)
 	    {
 	      if (u.ieee.exponent > v.ieee.exponent)
 		u.ieee.exponent += 2 * LDBL_MANT_DIG + 2;

-- 
Joseph S. Myers
joseph@codesourcery.com