[BZ #2423]

2006-03-07  Jakub Jelinek  <jakub@redhat.com>
	[BZ #2423]
	* math/libm-test.inc [TEST_LDOUBLE] (ceil_test, floor_test, rint_test,
	round_test, trunc_test): Only run some of the new tests if
	LDBL_MANT_DIG > 100.

2006-03-03  Steven Munroe  <sjmunroe@us.ibm.com>
	    Alan Modra  <amodra@bigpond.net.au>

	* sysdeps/powerpc/fpu/fenv_libc.h (__fegetround, __fesetround):
	Define inline implementations.
	* sysdeps/powerpc/fpu/fegetround.c: Use __fegetround.
	* sysdeps/powerpc/fpu/fesetround.c: Use __fesetround.

	* sysdeps/powerpc/fpu/math_ldbl.h: New file.

	[BZ #2423]
	* math/libm-test.inc [TEST_LDOUBLE] (ceil_test, floor_test, rint_test,
	round_test, trunc_test): Add new tests.
	* sysdeps/ieee754/ldbl-128ibm/math_ldbl.h
	(EXTRACT_IBM_EXTENDED_MANTISSA, INSERT_IBM_EXTENDED_MANTISSA):
	Removed, replaced with ...
	(ldbl_extract_mantissa, ldbl_insert_mantissa, ldbl_pack, ldbl_unpack,
	ldbl_canonicalise, ldbl_nearbyint): New functions.
	* sysdeps/ieee754/ldbl-128ibm/e_fmodl.c (__ieee754_fmodl): Replace
	EXTRACT_IBM_EXTENDED_MANTISSA and INSERT_IBM_EXTENDED_MANTISSA
	with ldbl_extract_mantissa and ldbl_insert_mantissa.
	* sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c (__ieee754_rem_pio2l):
	Replace EXTRACT_IBM_EXTENDED_MANTISSA with ldbl_extract_mantissa.
	(ldbl_extract_mantissa, ldbl_insert_mantissa): New inline functions.
	* sysdeps/ieee754/ldbl-128ibm/s_ceill.c (__ceill): Handle rounding
	that spans doubles in IBM long double format.
	* sysdeps/ieee754/ldbl-128ibm/s_floorl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/s_rintl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/s_roundl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/s_truncl.c: Likewise.
	* sysdeps/powerpc/powerpc64/fpu/s_rintl.S: File removed.

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);