From b9a76339be2514c700d801e179ef9b6c910eaedf Mon Sep 17 00:00:00 2001
From: Joseph Myers <joseph@codesourcery.com>
Date: Thu, 18 Feb 2016 22:24:32 +0000
Subject: [PATCH] Fix ldbl-128ibm roundl for non-default rounding modes (bug
 19594).

The ldbl-128ibm implementation of roundl is only correct in
round-to-nearest mode (in other modes, there are incorrect results and
overflow exceptions in some cases).  This patch reimplements it along
the lines used for floorl, ceill and truncl, using __round on the high
part, and on the low part if the high part is an integer, and then
adjusting in the cases where this is incorrect.

Tested for powerpc.

	[BZ #19594]
	* sysdeps/ieee754/ldbl-128ibm/s_roundl.c (__roundl): Use __round
	on high and low parts then adjust result and use
	ldbl_canonicalize_int if needed.
---
 ChangeLog                              |  5 ++
 sysdeps/ieee754/ldbl-128ibm/s_roundl.c | 70 +++++++++++++-------------
 2 files changed, 39 insertions(+), 36 deletions(-)

diff --git a/ChangeLog b/ChangeLog
index 1f7ba5cfae..f9c3118def 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,10 @@
 2016-02-18  Joseph Myers  <joseph@codesourcery.com>
 
+	[BZ #19594]
+	* sysdeps/ieee754/ldbl-128ibm/s_roundl.c (__roundl): Use __round
+	on high and low parts then adjust result and use
+	ldbl_canonicalize_int if needed.
+
 	[BZ #19593]
 	* sysdeps/ieee754/ldbl-128ibm/s_truncl.c (__truncl): Use __trunc
 	on high part and __floor or __ceil on low part then use
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_roundl.c b/sysdeps/ieee754/ldbl-128ibm/s_roundl.c
index 20813ed366..b01510fef3 100644
--- a/sysdeps/ieee754/ldbl-128ibm/s_roundl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_roundl.c
@@ -38,46 +38,44 @@ __roundl (long double x)
 			&& __builtin_isless (__builtin_fabs (xh),
 					     __builtin_inf ()), 1))
     {
-      double orig_xh;
-
-      /* Long double arithmetic, including the canonicalisation below,
-	 only works in round-to-nearest mode.  */
-
-      /* Convert the high double to integer.  */
-      orig_xh = xh;
-      hi = ldbl_nearbyint (xh);
-
-      /* Subtract integral high part from the value.  */
-      xh -= hi;
-      ldbl_canonicalize (&xh, &xl);
-
-      /* Now convert the low double, adjusted for any remainder from the
-	 high double.  */
-      lo = ldbl_nearbyint (xh);
-
-      /* Adjust the result when the remainder is exactly 0.5.  nearbyint
-	 rounds values halfway between integers to the nearest even
-	 integer.  roundl must round away from zero.
-	 Also correct cases where nearbyint returns an incorrect value
-	 for LO.  */
-      xh -= lo;
-      ldbl_canonicalize (&xh, &xl);
-      if (xh == 0.5)
+      hi = __round (xh);
+      if (hi != xh)
 	{
-	  if (xl > 0.0 || (xl == 0.0 && orig_xh > 0.0))
-	    lo += 1.0;
+	  /* The high part is not an integer; the low part only
+	     affects the result if the high part is exactly half way
+	     between two integers and the low part is nonzero with the
+	     opposite sign.  */
+	  if (fabs (hi - xh) == 0.5)
+	    {
+	      if (xh > 0 && xl < 0)
+		xh = hi - 1;
+	      else if (xh < 0 && xl > 0)
+		xh = hi + 1;
+	      else
+		xh = hi;
+	    }
+	  else
+	    xh = hi;
+	  xl = 0;
 	}
-      else if (-xh == 0.5)
+      else
 	{
-	  if (xl < 0.0 || (xl == 0.0 && orig_xh < 0.0))
-	    lo -= 1.0;
+	  /* The high part is a nonzero integer.  */
+	  lo = __round (xl);
+	  if (fabs (lo - xl) == 0.5)
+	    {
+	      if (xh > 0 && xl < 0)
+		xl = lo + 1;
+	      else if (xh < 0 && lo > 0)
+		xl = lo - 1;
+	      else
+		xl = lo;
+	    }
+	  else
+	    xl = lo;
+	  xh = hi;
+	  ldbl_canonicalize_int (&xh, &xl);
 	}
-
-      /* Ensure the final value is canonical.  In certain cases,
-	 rounding causes hi,lo calculated so far to be non-canonical.  */
-      xh = hi;
-      xl = lo;
-      ldbl_canonicalize (&xh, &xl);
     }
 
   return ldbl_pack (xh, xl);