Fix powerpc32 llround, llroundf exceptions (bug 19125).

The powerpc32 implementations of llroundf and llround produce spurious
and missing exceptions (some arising from such exceptions from
conversions to long long, some present even when fctidz is used).

This patch fixes those problems in a similar way to the llrint /
llrintf fixes.  The spurious exceptions in the fctidz case for large
arguments arise from a converted value that saturated as LLONG_MAX
being converted back to float or double (the conversion back being
inexact, but "inexact" must not be raised together with "invalid"),
and from the subtraction x - xrf also being inexact for sufficiently
large arguments (whether the saturation was to LLONG_MAX or
LLONG_MIN); those are fixed by returning early if the argument is
large enough that no rounding is needed.

This code is not used for --with-cpu=power4 builds (I suspect the code
used in that case may also produce spurious "inexact" exceptions, but
that's something to investigate later).

Tested for powerpc.

	[BZ #19125]
	* sysdeps/powerpc/powerpc32/fpu/s_llround.c: Include <limits.h>,
	<math_private.h> and <stdint.h>.
	(__llround): Avoid conversions to and from long long int, and
	subtractions, where those might raise spurious exceptions.
	* sysdeps/powerpc/powerpc32/fpu/s_llroundf.c: Include
	<math_private.h> and <stdint.h>.
	(__llroundf): Avoid conversions to and from long long int, and
	subtractions, where those might raise spurious exceptions.

2 files changed, 68 insertions, 2 deletions

diff --git a/sysdeps/powerpc/powerpc32/fpu/s_llround.c b/sysdeps/powerpc/powerpc32/fpu/s_llround.c
index 85514a1d3f..bdee74b13c 100644
--- a/sysdeps/powerpc/powerpc32/fpu/s_llround.c
+++ b/sysdeps/powerpc/powerpc32/fpu/s_llround.c
@@ -16,8 +16,11 @@
    License along with the GNU C Library; if not, see
    <http://www.gnu.org/licenses/>.  */
 
+#include <limits.h>
 #include <math.h>
 #include <math_ldbl_opt.h>
+#include <math_private.h>
+#include <stdint.h>
 
 /* Round to the nearest integer, with values exactly on a 0.5 boundary
    rounded away from zero, regardless of the current rounding mode.
@@ -27,7 +30,42 @@
 long long int
 __llround (double x)
 {
-  long long xr = (long long) x;
+  long long xr;
+  if (HAVE_PPC_FCTIDZ)
+    xr = (long long) x;
+  else
+    {
+      /* Avoid incorrect exceptions from libgcc conversions (as of GCC
+	 5): <https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59412>.  */
+      if (fabs (x) < 0x1p31)
+	xr = (long long int) (long int) x;
+      else
+	{
+	  uint64_t i0;
+	  EXTRACT_WORDS64 (i0, x);
+	  int exponent = ((i0 >> 52) & 0x7ff) - 0x3ff;
+	  if (exponent < 63)
+	    {
+	      unsigned long long int mant
+		= (i0 & ((1ULL << 52) - 1)) | (1ULL << 52);
+	      if (exponent < 52)
+		/* llround is not required to raise "inexact".  */
+		mant >>= 52 - exponent;
+	      else
+		mant <<= exponent - 52;
+	      xr = (long long int) ((i0 & (1ULL << 63)) != 0 ? -mant : mant);
+	    }
+	  else if (x == (double) LLONG_MIN)
+	    xr = LLONG_MIN;
+	  else
+	    xr = (long long int) (long int) x << 32;
+	}
+    }
+  /* Avoid spurious "inexact" converting LLONG_MAX to double, and from
+     subtraction when the result is out of range, by returning early
+     for arguments large enough that no rounding is needed.  */
+  if (!(fabs (x) < 0x1p52))
+    return xr;
   double xrf = (double) xr;
 
   if (x >= 0.0)
diff --git a/sysdeps/powerpc/powerpc32/fpu/s_llroundf.c b/sysdeps/powerpc/powerpc32/fpu/s_llroundf.c
index 8ce4b26191..6c67857ec5 100644
--- a/sysdeps/powerpc/powerpc32/fpu/s_llroundf.c
+++ b/sysdeps/powerpc/powerpc32/fpu/s_llroundf.c
@@ -17,6 +17,8 @@
    <http://www.gnu.org/licenses/>.  */
 
 #include <math.h>
+#include <math_private.h>
+#include <stdint.h>
 
 /* Round to the nearest integer, with values exactly on a 0.5 boundary
    rounded away from zero, regardless of the current rounding mode.
@@ -26,7 +28,33 @@
 long long int
 __llroundf (float x)
 {
-  long long xr = (long long) x;
+  long long xr;
+  if (HAVE_PPC_FCTIDZ)
+    xr = (long long) x;
+  else
+    {
+      float ax = fabsf (x);
+      /* Avoid incorrect exceptions from libgcc conversions (as of GCC
+	 5): <https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59412>.  */
+      if (ax < 0x1p31f)
+	xr = (long long int) (long int) x;
+      else if (!(ax < 0x1p55f))
+	xr = (long long int) (long int) (x * 0x1p-32f) << 32;
+      else
+	{
+	  uint32_t i0;
+	  GET_FLOAT_WORD (i0, x);
+	  int exponent = ((i0 >> 23) & 0xff) - 0x7f;
+	  unsigned long long int mant = (i0 & 0x7fffff) | 0x800000;
+	  mant <<= exponent - 23;
+	  xr = (long long int) ((i0 & 0x80000000) != 0 ? -mant : mant);
+	}
+    }
+  /* Avoid spurious "inexact" converting LLONG_MAX to float, and from
+     subtraction when the result is out of range, by returning early
+     for arguments large enough that no rounding is needed.  */
+  if (!(fabsf (x) < 0x1p23f))
+    return xr;
   float xrf = (float) xr;
 
   if (x >= 0.0)