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, 8 insertions, 9 deletions

diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index 34a63233d3..72a1f3c4c2 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -2069,15 +2069,14 @@ fma -min -min 0 missing-errno
 fma -min -min -0 missing-errno
 
 # Bug 6801: errno setting may be missing.
-# Bug 13304: results on directed rounding may be incorrect.
-fma max max min missing-errno xfail-rounding:ldbl-128ibm
-fma max max -min missing-errno xfail-rounding:ldbl-128ibm
-fma max -max min missing-errno xfail-rounding:ldbl-128ibm
-fma max -max -min missing-errno xfail-rounding:ldbl-128ibm
-fma -max max min missing-errno xfail-rounding:ldbl-128ibm
-fma -max max -min missing-errno xfail-rounding:ldbl-128ibm
-fma -max -max min missing-errno xfail-rounding:ldbl-128ibm
-fma -max -max -min missing-errno xfail-rounding:ldbl-128ibm
+fma max max min missing-errno
+fma max max -min missing-errno
+fma max -max min missing-errno
+fma max -max -min missing-errno
+fma -max max min missing-errno
+fma -max max -min missing-errno
+fma -max -max min missing-errno
+fma -max -max -min missing-errno
 
 fma 0x1.7ff8p+13 0x1.000002p+0 0x1.ffffp-24
 fma 0x1.fffp+0 0x1.00001p+0 -0x1.fffp+0