1 files changed, 51 insertions, 73 deletions

diff --git a/sysdeps/ieee754/ldbl-128ibm/s_roundl.c b/sysdeps/ieee754/ldbl-128ibm/s_roundl.c
index 7fe84b8..0880e6e 100644
--- a/sysdeps/ieee754/ldbl-128ibm/s_roundl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_roundl.c
@@ -22,7 +22,7 @@
    when it's coded in C.  */
 
 #include <math.h>
-#include <fenv.h>
+#include <fenv_libc.h>
 #include <math_ldbl_opt.h>
 #include <float.h>
 #include <ieee754.h>
@@ -37,84 +37,62 @@ __roundl (x)
      long double x;
 #endif
 {
-  static const double TWO52 = 4503599627370496.0;
-  static const double HALF = 0.5;
-  int mode = fegetround();
-  union ibm_extended_long_double u;
-  u.d = x;
-
-  if (fabs (u.dd[0]) < TWO52)
-    {      
-      fesetround(FE_TOWARDZERO);
-      if (u.dd[0] > 0.0)
-	{
-	  u.dd[0] += HALF;
-	  u.dd[0] += TWO52;
-	  u.dd[0] -= TWO52;
-	}
-      else if (u.dd[0] < 0.0)
-	{
-	  u.dd[0] = TWO52 - (u.dd[0] - HALF);
-	  u.dd[0] = -(u.dd[0] - TWO52);
-	}
-      u.dd[1] = 0.0;
-      fesetround(mode);
-    }
-  else if (fabs (u.dd[1]) < TWO52 && u.dd[1] != 0.0)
+  double xh, xl, hi, lo;
+
+  ldbl_unpack (x, &xh, &xl);
+
+  /* Return Inf, Nan, +/-0 unchanged.  */
+  if (__builtin_expect (xh != 0.0
+			&& __builtin_isless (__builtin_fabs (xh),
+					     __builtin_inf ()), 1))
     {
-      double high, low;
-      /* In this case we have to round the low double and handle any
-         adjustment to the high double that may be caused by rounding
-         (up).  This is complicated by the fact that the high double
-         may already be rounded and the low double may have the
-         opposite sign to compensate.  */
-      if (u.dd[0] > 0.0)
+      double orig_xh;
+      int save_round = fegetround ();
+
+      /* Long double arithmetic, including the canonicalisation below,
+	 only works in round-to-nearest mode.  */
+      fesetround (FE_TONEAREST);
+
+      /* 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)
 	{
-	  if (u.dd[1] > 0.0)
-	    {
-	      /* If the high/low doubles are the same sign then simply
-	         round the low double.  */
-	      high = u.dd[0];
-	      low = u.dd[1];
-	    }
-	  else if (u.dd[1] < 0.0)
-	    {
-	      /* Else the high double is pre rounded and we need to
-	         adjust for that.  */
-	      high = nextafter (u.dd[0], 0.0);
-	      low = u.dd[1] + (u.dd[0] - high);
-	    }     
-          fesetround(FE_TOWARDZERO);
-	  low += HALF;
-	  low += TWO52;
-	  low -= TWO52;
+	  if (xl > 0.0 || (xl == 0.0 && orig_xh > 0.0))
+	    lo += 1.0;
 	}
-      else if (u.dd[0] < 0.0)
+      else if (-xh == 0.5)
 	{
-	  if (u.dd[1] < 0.0)
-	    {
-	      /* If the high/low doubles are the same sign then simply
-	         round the low double.  */
-	      high = u.dd[0];
-	      low = u.dd[1];
-	    }
-	  else if (u.dd[1] > 0.0)
-	    {
-	      /* Else the high double is pre rounded and we need to
-	         adjust for that.  */
-	      high = nextafter (u.dd[0], 0.0);
-	      low = u.dd[1] + (u.dd[0] - high);
-	    }     
-          fesetround(FE_TOWARDZERO);
-	  low -= HALF;
-	  low = TWO52 - low;
-	  low = -(low - TWO52);
+	  if (xl < 0.0 || (xl == 0.0 && orig_xh < 0.0))
+	    lo -= 1.0;
 	}
-      fesetround(mode);
-      u.dd[0] = high + low;
-      u.dd[1] = high - u.dd[0] + low;
+
+      /* 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);
+
+      fesetround (save_round);
     }
-  return u.d;
+
+  return ldbl_pack (xh, xl);
 }
 
 long_double_symbol (libm, __roundl, roundl);