Implement proper fmal for ldbl-128ibm (bug 13304).

ldbl-128ibm had an implementation of fmal that just did (x * y) + z in
most cases, with no attempt at actually being a fused operation.

This patch replaces it with a genuine fused operation.  It is not
necessarily correctly rounding, but should produce a result at least
as accurate as the long double arithmetic operations in libgcc, which
I think is all that can reasonably be expected for such a non-IEEE
format where arithmetic is approximate rather than rounded according
to any particular rule for determining the exact result.  Like the
libgcc arithmetic, it may produce spurious overflow and underflow
results, and it falls back to the libgcc multiplication in the case of
(finite, finite, zero).

This concludes the fixes for bug 13304; any subsequently found fma
issues should go in separate Bugzilla bugs.  Various other pieces of
bug 13304 were fixed in past releases over the past several years.

Tested for powerpc.

	[BZ #13304]
	* sysdeps/ieee754/ldbl-128ibm/s_fmal.c: Include <fenv.h>,
	<float.h>, <math_private.h> and <stdlib.h>.
	(add_split): New function.
	(mul_split): Likewise.
	(ext_val): New typedef.
	(store_ext_val): New function.
	(mul_ext_val): New function.
	(compare): New function.
	(add_split_ext): New function.
	(__fmal): After checking for Inf, NaN and zero, compute result as
	an exact sum of scaled double values in round-to-nearest before
	adding those up and adjusting for other rounding modes.
	* math/auto-libm-test-in: Remove xfail-rounding:ldbl-128ibm from
	tests of fma.
	* math/auto-libm-test-out: Regenerated.

1 files changed, 249 insertions, 11 deletions

diff --git a/sysdeps/ieee754/ldbl-128ibm/s_fmal.c b/sysdeps/ieee754/ldbl-128ibm/s_fmal.c
index eb3ee3cfb8..177a04817b 100644
--- a/sysdeps/ieee754/ldbl-128ibm/s_fmal.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_fmal.c
@@ -17,25 +17,263 @@
    License along with the GNU C Library; if not, see
    <http://www.gnu.org/licenses/>.  */
 
+#include <fenv.h>
+#include <float.h>
 #include <math.h>
+#include <math_private.h>
 #include <math_ldbl_opt.h>
+#include <stdlib.h>
+
+/* Calculate X + Y exactly and store the result in *HI + *LO.  It is
+   given that |X| >= |Y| and the values are small enough that no
+   overflow occurs.  */
+
+static void
+add_split (double *hi, double *lo, double x, double y)
+{
+  /* Apply Dekker's algorithm.  */
+  *hi = x + y;
+  *lo = (x - *hi) + y;
+}
+
+/* Calculate X * Y exactly and store the result in *HI + *LO.  It is
+   given that the values are small enough that no overflow occurs and
+   large enough (or zero) that no underflow occurs.  */
+
+static void
+mul_split (double *hi, double *lo, double x, double y)
+{
+#ifdef __FP_FAST_FMA
+  /* Fast built-in fused multiply-add.  */
+  *hi = x * y;
+  *lo = __builtin_fma (x, y, -*hi);
+#else
+  /* Apply Dekker's algorithm.  */
+  *hi = x * y;
+# define C ((1 << (DBL_MANT_DIG + 1) / 2) + 1)
+  double x1 = x * C;
+  double y1 = y * C;
+# undef C
+  x1 = (x - x1) + x1;
+  y1 = (y - y1) + y1;
+  double x2 = x - x1;
+  double y2 = y - y1;
+  *lo = (((x1 * y1 - *hi) + x1 * y2) + x2 * y1) + x2 * y2;
+#endif
+}
+
+/* Value with extended range, used in intermediate computations.  */
+typedef struct
+{
+  /* Value in [0.5, 1), as from frexp, or 0.  */
+  double val;
+  /* Exponent of power of 2 it is multiplied by, or 0 for zero.  */
+  int exp;
+} ext_val;
+
+/* Store D as an ext_val value.  */
+
+static void
+store_ext_val (ext_val *v, double d)
+{
+  v->val = __frexp (d, &v->exp);
+}
+
+/* Store X * Y as ext_val values *V0 and *V1.  */
+
+static void
+mul_ext_val (ext_val *v0, ext_val *v1, double x, double y)
+{
+  int xexp, yexp;
+  x = __frexp (x, &xexp);
+  y = __frexp (y, &yexp);
+  double hi, lo;
+  mul_split (&hi, &lo, x, y);
+  store_ext_val (v0, hi);
+  if (hi != 0)
+    v0->exp += xexp + yexp;
+  store_ext_val (v1, lo);
+  if (lo != 0)
+    v1->exp += xexp + yexp;
+}
+
+/* Compare absolute values of ext_val values pointed to by P and Q for
+   qsort.  */
+
+static int
+compare (const void *p, const void *q)
+{
+  const ext_val *pe = p;
+  const ext_val *qe = q;
+  if (pe->val == 0)
+    return qe->val == 0 ? 0 : -1;
+  else if (qe->val == 0)
+    return 1;
+  else if (pe->exp < qe->exp)
+    return -1;
+  else if (pe->exp > qe->exp)
+    return 1;
+  else
+    {
+      double pd = fabs (pe->val);
+      double qd = fabs (qe->val);
+      if (pd < qd)
+	return -1;
+      else if (pd == qd)
+	return 0;
+      else
+	return 1;
+    }
+}
+
+/* Calculate *X + *Y exactly, storing the high part in *X (rounded to
+   nearest) and the low part in *Y.  It is given that |X| >= |Y|.  */
+
+static void
+add_split_ext (ext_val *x, ext_val *y)
+{
+  int xexp = x->exp, yexp = y->exp;
+  if (y->val == 0 || xexp - yexp > 53)
+    return;
+  double hi = x->val;
+  double lo = __scalbn (y->val, yexp - xexp);
+  add_split (&hi, &lo, hi, lo);
+  store_ext_val (x, hi);
+  if (hi != 0)
+    x->exp += xexp;
+  store_ext_val (y, lo);
+  if (lo != 0)
+    y->exp += xexp;
+}
 
 long double
 __fmal (long double x, long double y, long double z)
 {
-	/* An IBM long double 128 is really just 2 IEEE64 doubles, and in
-	 * the case of inf/nan only the first double counts. So we use the
-	 * (double) cast to avoid any data movement.   */
-       if ((isfinite ((double)x) && isfinite ((double)y)) && isinf ((double)z))
-               return (z);
+  double xhi, xlo, yhi, ylo, zhi, zlo;
+  int64_t hx, hy, hz;
+  int xexp, yexp, zexp;
+  double scale_val;
+  int scale_exp;
+  ldbl_unpack (x, &xhi, &xlo);
+  EXTRACT_WORDS64 (hx, xhi);
+  xexp = (hx & 0x7ff0000000000000LL) >> 52;
+  ldbl_unpack (y, &yhi, &ylo);
+  EXTRACT_WORDS64 (hy, yhi);
+  yexp = (hy & 0x7ff0000000000000LL) >> 52;
+  ldbl_unpack (z, &zhi, &zlo);
+  EXTRACT_WORDS64 (hz, zhi);
+  zexp = (hz & 0x7ff0000000000000LL) >> 52;
+
+  /* If z is Inf or NaN, but x and y are finite, avoid any exceptions
+     from computing x * y.  */
+  if (zexp == 0x7ff && xexp != 0x7ff && yexp != 0x7ff)
+    return (z + x) + y;
+
+  /* If z is zero and x are y are nonzero, compute the result as x * y
+     to avoid the wrong sign of a zero result if x * y underflows to
+     0.  */
+  if (z == 0 && x != 0 && y != 0)
+    return x * y;
+
+  /* If x or y or z is Inf/NaN, or if x * y is zero, compute as x * y
+     + z.  */
+  if (xexp == 0x7ff || yexp == 0x7ff || zexp == 0x7ff
+      || x == 0 || y == 0)
+    return (x * y) + z;
+
+  {
+    SET_RESTORE_ROUND (FE_TONEAREST);
+
+    ext_val vals[10];
+    store_ext_val (&vals[0], zhi);
+    store_ext_val (&vals[1], zlo);
+    mul_ext_val (&vals[2], &vals[3], xhi, yhi);
+    mul_ext_val (&vals[4], &vals[5], xhi, ylo);
+    mul_ext_val (&vals[6], &vals[7], xlo, yhi);
+    mul_ext_val (&vals[8], &vals[9], xlo, ylo);
+    qsort (vals, 10, sizeof (ext_val), compare);
+    /* Add up the values so that each element of VALS has absolute
+       value at most equal to the last set bit of the next nonzero
+       element.  */
+    for (size_t i = 0; i <= 8; i++)
+      {
+	add_split_ext (&vals[i + 1], &vals[i]);
+	qsort (vals + i + 1, 9 - i, sizeof (ext_val), compare);
+      }
+    /* Add up the values in the other direction, so that each element
+       of VALS has absolute value less than 5ulp of the next
+       value.  */
+    size_t dstpos = 9;
+    for (size_t i = 1; i <= 9; i++)
+      {
+	if (vals[dstpos].val == 0)
+	  {
+	    vals[dstpos] = vals[9 - i];
+	    vals[9 - i].val = 0;
+	    vals[9 - i].exp = 0;
+	  }
+	else
+	  {
+	    add_split_ext (&vals[dstpos], &vals[9 - i]);
+	    if (vals[9 - i].val != 0)
+	      {
+		if (9 - i < dstpos - 1)
+		  {
+		    vals[dstpos - 1] = vals[9 - i];
+		    vals[9 - i].val = 0;
+		    vals[9 - i].exp = 0;
+		  }
+		dstpos--;
+	      }
+	  }
+      }
+    /* If the result is an exact zero, it results from adding two
+       values with opposite signs; recompute in the original rounding
+       mode.  */
+    if (vals[9].val == 0)
+      goto zero_out;
+    /* Adding the top three values will now give a result as accurate
+       as the underlying long double arithmetic.  */
+    add_split_ext (&vals[9], &vals[8]);
+    if (compare (&vals[8], &vals[7]) < 0)
+      {
+	ext_val tmp = vals[7];
+	vals[7] = vals[8];
+	vals[8] = tmp;
+      }
+    add_split_ext (&vals[8], &vals[7]);
+    add_split_ext (&vals[9], &vals[8]);
+    if (vals[9].exp > DBL_MAX_EXP || vals[9].exp < DBL_MIN_EXP)
+      {
+	/* Overflow or underflow, with the result depending on the
+	   original rounding mode, but not on the low part computed
+	   here.  */
+	scale_val = vals[9].val;
+	scale_exp = vals[9].exp;
+	goto scale_out;
+      }
+    double hi = __scalbn (vals[9].val, vals[9].exp);
+    double lo = __scalbn (vals[8].val, vals[8].exp);
+    /* It is possible that the low part became subnormal and was
+       rounded so that the result is no longer canonical.  */
+    ldbl_canonicalize (&hi, &lo);
+    long double ret = ldbl_pack (hi, lo);
+    math_check_force_underflow (ret);
+    return ret;
+  }
 
-       /* If z is zero and x are y are nonzero, compute the result
-	  as x * y to avoid the wrong sign of a zero result if x * y
-	  underflows to 0.  */
-       if (z == 0 && x != 0 && y != 0)
-	 return x * y;
+ scale_out:
+  scale_val = math_opt_barrier (scale_val);
+  scale_val = __scalbn (scale_val, scale_exp);
+  if (fabs (scale_val) == DBL_MAX)
+    return __copysignl (LDBL_MAX, scale_val);
+  math_check_force_underflow (scale_val);
+  return scale_val;
 
-       return (x * y) + z;
+ zero_out:;
+  double zero = 0.0;
+  zero = math_opt_barrier (zero);
+  return zero - zero;
 }
 #if IS_IN (libm)
 long_double_symbol (libm, __fmal, fmal);