This is the mail archive of the libc-alpha@sourceware.org mailing list for the glibc project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

[PATCH] Fix long double rounding functions


With additional testing we found more corner cases where the ldbl-128ibm
implementations of the rounding functions gave incorrect results. For
example:

Test "rint (4503599627370494.5000000000001) == 4503599627370495.0":
ldouble: 18014398509481984
Failure: Test: rint (4503599627370494.5000000000001) == 4503599627370495.0
Result:
 is:          4.50359962737049400000e+15  
0x1.ffffffffffffc00000000000000p+51
 should be:   4.50359962737049500000e+15  
0x1.ffffffffffffe00000000000000p+51
 difference:  1.00000000000000000000e+00  
0x1.000000000000000000000000000p+0
 ulp       :  18014398509481984.0000
 max.ulp   :  0.0000
Test "rint (-4503599627370494.5000000000001) == -4503599627370495.0":
ldouble: 18014398509481984
Failure: Test: rint (-4503599627370494.5000000000001) == -4503599627370495.0
Result:
 is:         -4.50359962737049400000e+15 
-0x1.ffffffffffffc00000000000000p+51
 should be:  -4.50359962737049500000e+15 
-0x1.ffffffffffffe00000000000000p+51
 difference:  1.00000000000000000000e+00  
0x1.000000000000000000000000000p+0
 ulp       :  18014398509481984.0000
 max.ulp   :  0.0000

These are cases where the "decimal point" is at or near the split
between the high/low doubles that compose the long double. Working with
Alan Modra we developed more comprehensive tests and improved algorithms
for rintl, roundl, ceill, floorl, truncl.

These improved implementation exposed an performance issue in for
powerpc due to  the heavy use of fegetround/fesetround. So we provided
an improved implementation of fesetround and provided for inlining
fegetround/fesetround for the long double rounding functions.

Finally we provided updated tests for libm-test.inc.

The remaining work is to improve llrintl and llroundl which we will
provide soon as a separate patch.

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

	* math/libm-test.inc [TEST_LDOUBLE] (ceil_test, floor_test, rint_test,
	round_test, trunc_test): Add new tests.
	
	* sysdeps/powerpc/fpu/fenv_libc.h (___fegetround, ___fesetround):
	Define inline implementations.
	* sysdeps/powerpc/fpu/fesetround.c: Improved implementation.
	
	* sysdeps/ieee754/ldbl-128ibm/math_ldbl.h (__ldbl_pack, __ldbl_unpack,
	__ldbl_canonicalise, __ldbl_nearbyint): Define inline utility
	functions for IBM long double format.
	* 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/fpu/math_ldbl.h: New file.
	* sysdeps/powerpc/powerpc64/fpu/s_rintl.S: Removed.

diff -urN libc24-cvstip-20060301/math/libm-test.inc libc24/math/libm-test.inc
--- libc24-cvstip-20060301/math/libm-test.inc	2006-01-27 18:35:55.000000000 -0600
+++ libc24/math/libm-test.inc	2006-03-01 17:05:50.000000000 -0600
@@ -1628,8 +1628,12 @@
 
   TEST_f_f (ceil, M_PIl, 4.0);
   TEST_f_f (ceil, -M_PIl, -3.0);
+  TEST_f_f (ceil, 0.1, 1.0);
   TEST_f_f (ceil, 0.25, 1.0);
+  TEST_f_f (ceil, 0.625, 1.0);
+  TEST_f_f (ceil, -0.1, minus_zero);
   TEST_f_f (ceil, -0.25, minus_zero);
+  TEST_f_f (ceil, -0.625, minus_zero);
 
 #ifdef TEST_LDOUBLE
   /* The result can only be represented in long double.  */
@@ -1644,6 +1648,13 @@
   TEST_f_f (ceil, -4503599627370496.5L, -4503599627370496.0L);
   TEST_f_f (ceil, -4503599627370496.75L, -4503599627370496.0L);
   TEST_f_f (ceil, -4503599627370497.5L, -4503599627370497.0L);
+  
+  TEST_f_f (ceil, 4503599627370494.5000000000001L, 4503599627370495.0L);
+  TEST_f_f (ceil, 4503599627370495.5000000000001L, 4503599627370496.0L);
+  TEST_f_f (ceil, 4503599627370496.5000000000001L, 4503599627370497.0L);
+  TEST_f_f (ceil, -4503599627370494.5000000000001L, -4503599627370494.0L);
+  TEST_f_f (ceil, -4503599627370495.5000000000001L, -4503599627370495.0L);
+  TEST_f_f (ceil, -4503599627370496.5000000000001L, -4503599627370496.0L);
 
   TEST_f_f (ceil, 9007199254740991.5L, 9007199254740992.0L);
   TEST_f_f (ceil, 9007199254740992.25L, 9007199254740993.0L);
@@ -1657,6 +1668,20 @@
   TEST_f_f (ceil, -9007199254740992.75L, -9007199254740992.0L);
   TEST_f_f (ceil, -9007199254740993.5L, -9007199254740993.0L);
 
+  TEST_f_f (ceil, 9007199254740991.0000000000001L, 9007199254740992.0L);
+  TEST_f_f (ceil, 9007199254740992.0000000000001L, 9007199254740993.0L);
+  TEST_f_f (ceil, 9007199254740993.0000000000001L, 9007199254740994.0L);
+  TEST_f_f (ceil, 9007199254740991.5000000000001L, 9007199254740992.0L);
+  TEST_f_f (ceil, 9007199254740992.5000000000001L, 9007199254740993.0L);
+  TEST_f_f (ceil, 9007199254740993.5000000000001L, 9007199254740994.0L);
+
+  TEST_f_f (ceil, -9007199254740991.0000000000001L, -9007199254740991.0L);
+  TEST_f_f (ceil, -9007199254740992.0000000000001L, -9007199254740992.0L);
+  TEST_f_f (ceil, -9007199254740993.0000000000001L, -9007199254740993.0L);
+  TEST_f_f (ceil, -9007199254740991.5000000000001L, -9007199254740991.0L);
+  TEST_f_f (ceil, -9007199254740992.5000000000001L, -9007199254740992.0L);
+  TEST_f_f (ceil, -9007199254740993.5000000000001L, -9007199254740993.0L);
+
   TEST_f_f (ceil, 72057594037927935.5L, 72057594037927936.0L);
   TEST_f_f (ceil, 72057594037927936.25L, 72057594037927937.0L);
   TEST_f_f (ceil, 72057594037927936.5L, 72057594037927937.0L);
@@ -2628,9 +2653,12 @@
   TEST_f_f (floor, M_PIl, 3.0);
   TEST_f_f (floor, -M_PIl, -4.0);
 
+  TEST_f_f (floor, 0.1, 0.0);
   TEST_f_f (floor, 0.25, 0.0);
+  TEST_f_f (floor, 0.625, 0.0);
+  TEST_f_f (floor, -0.1, -1.0);
   TEST_f_f (floor, -0.25, -1.0);
-
+  TEST_f_f (floor, -0.625, -1.0);
 
 #ifdef TEST_LDOUBLE
   /* The result can only be represented in long double.  */
@@ -2639,12 +2667,18 @@
   TEST_f_f (floor, 4503599627370496.5L, 4503599627370496.0L);
   TEST_f_f (floor, 4503599627370496.75L, 4503599627370496.0L);
   TEST_f_f (floor, 4503599627370497.5L, 4503599627370497.0L);
+  TEST_f_f (floor, 4503599627370494.5000000000001L, 4503599627370494.0L);
+  TEST_f_f (floor, 4503599627370495.5000000000001L, 4503599627370495.0L);
+  TEST_f_f (floor, 4503599627370496.5000000000001L, 4503599627370496.0L);
 
   TEST_f_f (floor, -4503599627370495.5L, -4503599627370496.0L);
   TEST_f_f (floor, -4503599627370496.25L, -4503599627370497.0L);
   TEST_f_f (floor, -4503599627370496.5L, -4503599627370497.0L);
   TEST_f_f (floor, -4503599627370496.75L, -4503599627370497.0L);
   TEST_f_f (floor, -4503599627370497.5L, -4503599627370498.0L);
+  TEST_f_f (floor, -4503599627370494.5000000000001L, -4503599627370495.0L);
+  TEST_f_f (floor, -4503599627370495.5000000000001L, -4503599627370496.0L);
+  TEST_f_f (floor, -4503599627370496.5000000000001L, -4503599627370497.0L);
 
   TEST_f_f (floor, 9007199254740991.5L, 9007199254740991.0L);
   TEST_f_f (floor, 9007199254740992.25L, 9007199254740992.0L);
@@ -2652,12 +2686,26 @@
   TEST_f_f (floor, 9007199254740992.75L, 9007199254740992.0L);
   TEST_f_f (floor, 9007199254740993.5L, 9007199254740993.0L);
 
+  TEST_f_f (floor, 9007199254740991.0000000000001L, 9007199254740991.0L);
+  TEST_f_f (floor, 9007199254740992.0000000000001L, 9007199254740992.0L);
+  TEST_f_f (floor, 9007199254740993.0000000000001L, 9007199254740993.0L);
+  TEST_f_f (floor, 9007199254740991.5000000000001L, 9007199254740991.0L);
+  TEST_f_f (floor, 9007199254740992.5000000000001L, 9007199254740992.0L);
+  TEST_f_f (floor, 9007199254740993.5000000000001L, 9007199254740993.0L);
+
   TEST_f_f (floor, -9007199254740991.5L, -9007199254740992.0L);
   TEST_f_f (floor, -9007199254740992.25L, -9007199254740993.0L);
   TEST_f_f (floor, -9007199254740992.5L, -9007199254740993.0L);
   TEST_f_f (floor, -9007199254740992.75L, -9007199254740993.0L);
   TEST_f_f (floor, -9007199254740993.5L, -9007199254740994.0L);
 
+  TEST_f_f (floor, -9007199254740991.0000000000001L, -9007199254740992.0L);
+  TEST_f_f (floor, -9007199254740992.0000000000001L, -9007199254740993.0L);
+  TEST_f_f (floor, -9007199254740993.0000000000001L, -9007199254740994.0L);
+  TEST_f_f (floor, -9007199254740991.5000000000001L, -9007199254740992.0L);
+  TEST_f_f (floor, -9007199254740992.5000000000001L, -9007199254740993.0L);
+  TEST_f_f (floor, -9007199254740993.5000000000001L, -9007199254740994.0L);
+
   TEST_f_f (floor, 72057594037927935.5L, 72057594037927935.0L);
   TEST_f_f (floor, 72057594037927936.25L, 72057594037927936.0L);
   TEST_f_f (floor, 72057594037927936.5L, 72057594037927936.0L);
@@ -3971,6 +4019,12 @@
   TEST_f_f (rint, -2.5, -2.0);
   TEST_f_f (rint, -3.5, -4.0);
   TEST_f_f (rint, -4.5, -4.0);
+  TEST_f_f (rint, 0.1, 0.0);
+  TEST_f_f (rint, 0.25, 0.0);
+  TEST_f_f (rint, 0.625, 1.0);
+  TEST_f_f (rint, -0.1, -0.0);
+  TEST_f_f (rint, -0.25, -0.0);
+  TEST_f_f (rint, -0.625, -1.0);
 #ifdef TEST_LDOUBLE
   /* The result can only be represented in long double.  */
   TEST_f_f (rint, 4503599627370495.5L, 4503599627370496.0L);
@@ -3978,12 +4032,34 @@
   TEST_f_f (rint, 4503599627370496.5L, 4503599627370496.0L);
   TEST_f_f (rint, 4503599627370496.75L, 4503599627370497.0L);
   TEST_f_f (rint, 4503599627370497.5L, 4503599627370498.0L);
+  
+  TEST_f_f (rint, 4503599627370494.5000000000001L, 4503599627370495.0L);
+  TEST_f_f (rint, 4503599627370495.5000000000001L, 4503599627370496.0L);
+  TEST_f_f (rint, 4503599627370496.5000000000001L, 4503599627370497.0L);
 
   TEST_f_f (rint, -4503599627370495.5L, -4503599627370496.0L);
   TEST_f_f (rint, -4503599627370496.25L, -4503599627370496.0L);
   TEST_f_f (rint, -4503599627370496.5L, -4503599627370496.0L);
   TEST_f_f (rint, -4503599627370496.75L, -4503599627370497.0L);
   TEST_f_f (rint, -4503599627370497.5L, -4503599627370498.0L);
+  
+  TEST_f_f (rint, -4503599627370494.5000000000001L, -4503599627370495.0L);
+  TEST_f_f (rint, -4503599627370495.5000000000001L, -4503599627370496.0L);
+  TEST_f_f (rint, -4503599627370496.5000000000001L, -4503599627370497.0L);
+
+  TEST_f_f (rint, 9007199254740991.0000000000001L, 9007199254740991.0L);
+  TEST_f_f (rint, 9007199254740992.0000000000001L, 9007199254740992.0L);
+  TEST_f_f (rint, 9007199254740993.0000000000001L, 9007199254740993.0L);
+  TEST_f_f (rint, 9007199254740991.5000000000001L, 9007199254740992.0L);
+  TEST_f_f (rint, 9007199254740992.5000000000001L, 9007199254740993.0L);
+  TEST_f_f (rint, 9007199254740993.5000000000001L, 9007199254740994.0L);
+
+  TEST_f_f (rint, -9007199254740991.0000000000001L, -9007199254740991.0L);
+  TEST_f_f (rint, -9007199254740992.0000000000001L, -9007199254740992.0L);
+  TEST_f_f (rint, -9007199254740993.0000000000001L, -9007199254740993.0L);
+  TEST_f_f (rint, -9007199254740991.5000000000001L, -9007199254740992.0L);
+  TEST_f_f (rint, -9007199254740992.5000000000001L, -9007199254740993.0L);
+  TEST_f_f (rint, -9007199254740993.5000000000001L, -9007199254740994.0L);
 
   TEST_f_f (rint, 9007199254740991.5L, 9007199254740992.0L);
   TEST_f_f (rint, 9007199254740992.25L, 9007199254740992.0L);
@@ -4039,6 +4115,45 @@
     TEST_f_f (rint, -1.0, -1.0);
     TEST_f_f (rint, -1.5, -2.0);
     TEST_f_f (rint, -2.0, -2.0);
+    TEST_f_f (rint, 0.1, 0.0);
+    TEST_f_f (rint, 0.25, 0.0);
+    TEST_f_f (rint, 0.625, 1.0);
+    TEST_f_f (rint, -0.1, -0.0);
+    TEST_f_f (rint, -0.25, -0.0);
+    TEST_f_f (rint, -0.625, -1.0);
+#ifdef TEST_LDOUBLE
+  /* The result can only be represented in long double.  */
+    TEST_f_f (rint, 4503599627370495.5L, 4503599627370496.0L);
+    TEST_f_f (rint, 4503599627370496.25L, 4503599627370496.0L);
+    TEST_f_f (rint, 4503599627370496.5L, 4503599627370496.0L);
+    TEST_f_f (rint, 4503599627370496.75L, 4503599627370497.0L);
+    TEST_f_f (rint, 4503599627370497.5L, 4503599627370498.0L);
+    TEST_f_f (rint, 4503599627370494.5000000000001L, 4503599627370495.0L);
+    TEST_f_f (rint, 4503599627370495.5000000000001L, 4503599627370496.0L);
+    TEST_f_f (rint, 4503599627370496.5000000000001L, 4503599627370497.0L);
+    TEST_f_f (rint, -4503599627370495.5L, -4503599627370496.0L);
+    TEST_f_f (rint, -4503599627370496.25L, -4503599627370496.0L);
+    TEST_f_f (rint, -4503599627370496.5L, -4503599627370496.0L);
+    TEST_f_f (rint, -4503599627370496.75L, -4503599627370497.0L);
+    TEST_f_f (rint, -4503599627370497.5L, -4503599627370498.0L);
+    TEST_f_f (rint, -4503599627370494.5000000000001L, -4503599627370495.0L);
+    TEST_f_f (rint, -4503599627370495.5000000000001L, -4503599627370496.0L);
+    TEST_f_f (rint, -4503599627370496.5000000000001L, -4503599627370497.0L);
+
+    TEST_f_f (rint, 9007199254740991.0000000000001L, 9007199254740991.0L);
+    TEST_f_f (rint, 9007199254740992.0000000000001L, 9007199254740992.0L);
+    TEST_f_f (rint, 9007199254740993.0000000000001L, 9007199254740993.0L);
+    TEST_f_f (rint, 9007199254740991.5000000000001L, 9007199254740992.0L);
+    TEST_f_f (rint, 9007199254740992.5000000000001L, 9007199254740993.0L);
+    TEST_f_f (rint, 9007199254740993.5000000000001L, 9007199254740994.0L);
+
+    TEST_f_f (rint, -9007199254740991.0000000000001L, -9007199254740991.0L);
+    TEST_f_f (rint, -9007199254740992.0000000000001L, -9007199254740992.0L);
+    TEST_f_f (rint, -9007199254740993.0000000000001L, -9007199254740993.0L);
+    TEST_f_f (rint, -9007199254740991.5000000000001L, -9007199254740992.0L);
+    TEST_f_f (rint, -9007199254740992.5000000000001L, -9007199254740993.0L);
+    TEST_f_f (rint, -9007199254740993.5000000000001L, -9007199254740994.0L);
+#endif
   }
 
   fesetround(save_round_mode);
@@ -4066,6 +4181,45 @@
     TEST_f_f (rint, -1.0, -1.0);
     TEST_f_f (rint, -1.5, -1.0);
     TEST_f_f (rint, -2.0, -2.0);
+    TEST_f_f (rint, 0.1, 0.0);
+    TEST_f_f (rint, 0.25, 0.0);
+    TEST_f_f (rint, 0.625, 0.0);
+    TEST_f_f (rint, -0.1, -0.0);
+    TEST_f_f (rint, -0.25, -0.0);
+    TEST_f_f (rint, -0.625, -0.0);
+#ifdef TEST_LDOUBLE
+  /* The result can only be represented in long double.  */
+    TEST_f_f (rint, 4503599627370495.5L, 4503599627370495.0L);
+    TEST_f_f (rint, 4503599627370496.25L, 4503599627370496.0L);
+    TEST_f_f (rint, 4503599627370496.5L, 4503599627370496.0L);
+    TEST_f_f (rint, 4503599627370496.75L, 4503599627370496.0L);
+    TEST_f_f (rint, 4503599627370497.5L, 4503599627370497.0L);
+    TEST_f_f (rint, 4503599627370494.5000000000001L, 4503599627370494.0L);
+    TEST_f_f (rint, 4503599627370495.5000000000001L, 4503599627370495.0L);
+    TEST_f_f (rint, 4503599627370496.5000000000001L, 4503599627370496.0L);
+    TEST_f_f (rint, -4503599627370495.5L, -4503599627370495.0L);
+    TEST_f_f (rint, -4503599627370496.25L, -4503599627370496.0L);
+    TEST_f_f (rint, -4503599627370496.5L, -4503599627370496.0L);
+    TEST_f_f (rint, -4503599627370496.75L, -4503599627370496.0L);
+    TEST_f_f (rint, -4503599627370497.5L, -4503599627370497.0L);
+    TEST_f_f (rint, -4503599627370494.5000000000001L, -4503599627370494.0L);
+    TEST_f_f (rint, -4503599627370495.5000000000001L, -4503599627370495.0L);
+    TEST_f_f (rint, -4503599627370496.5000000000001L, -4503599627370496.0L);
+
+    TEST_f_f (rint, 9007199254740991.0000000000001L, 9007199254740991.0L);
+    TEST_f_f (rint, 9007199254740992.0000000000001L, 9007199254740992.0L);
+    TEST_f_f (rint, 9007199254740993.0000000000001L, 9007199254740993.0L);
+    TEST_f_f (rint, 9007199254740991.5000000000001L, 9007199254740991.0L);
+    TEST_f_f (rint, 9007199254740992.5000000000001L, 9007199254740992.0L);
+    TEST_f_f (rint, 9007199254740993.5000000000001L, 9007199254740993.0L);
+
+    TEST_f_f (rint, -9007199254740991.0000000000001L, -9007199254740991.0L);
+    TEST_f_f (rint, -9007199254740992.0000000000001L, -9007199254740992.0L);
+    TEST_f_f (rint, -9007199254740993.0000000000001L, -9007199254740993.0L);
+    TEST_f_f (rint, -9007199254740991.5000000000001L, -9007199254740991.0L);
+    TEST_f_f (rint, -9007199254740992.5000000000001L, -9007199254740992.0L);
+    TEST_f_f (rint, -9007199254740993.5000000000001L, -9007199254740993.0L);
+#endif
   }
 
   fesetround(save_round_mode);
@@ -4093,6 +4247,45 @@
     TEST_f_f (rint, -1.0, -1.0);
     TEST_f_f (rint, -1.5, -2.0);
     TEST_f_f (rint, -2.0, -2.0);
+    TEST_f_f (rint, 0.1, 0.0);
+    TEST_f_f (rint, 0.25, 0.0);
+    TEST_f_f (rint, 0.625, 0.0);
+    TEST_f_f (rint, -0.1, -1.0);
+    TEST_f_f (rint, -0.25, -1.0);
+    TEST_f_f (rint, -0.625, -1.0);
+#ifdef TEST_LDOUBLE
+  /* The result can only be represented in long double.  */
+    TEST_f_f (rint, 4503599627370495.5L, 4503599627370495.0L);
+    TEST_f_f (rint, 4503599627370496.25L, 4503599627370496.0L);
+    TEST_f_f (rint, 4503599627370496.5L, 4503599627370496.0L);
+    TEST_f_f (rint, 4503599627370496.75L, 4503599627370496.0L);
+    TEST_f_f (rint, 4503599627370497.5L, 4503599627370497.0L);
+    TEST_f_f (rint, 4503599627370494.5000000000001L, 4503599627370494.0L);
+    TEST_f_f (rint, 4503599627370495.5000000000001L, 4503599627370495.0L);
+    TEST_f_f (rint, 4503599627370496.5000000000001L, 4503599627370496.0L);
+    TEST_f_f (rint, -4503599627370495.5L, -4503599627370496.0L);
+    TEST_f_f (rint, -4503599627370496.25L, -4503599627370497.0L);
+    TEST_f_f (rint, -4503599627370496.5L, -4503599627370497.0L);
+    TEST_f_f (rint, -4503599627370496.75L, -4503599627370497.0L);
+    TEST_f_f (rint, -4503599627370497.5L, -4503599627370498.0L);
+    TEST_f_f (rint, -4503599627370494.5000000000001L, -4503599627370495.0L);
+    TEST_f_f (rint, -4503599627370495.5000000000001L, -4503599627370496.0L);
+    TEST_f_f (rint, -4503599627370496.5000000000001L, -4503599627370497.0L);
+
+    TEST_f_f (rint, 9007199254740991.0000000000001L, 9007199254740991.0L);
+    TEST_f_f (rint, 9007199254740992.0000000000001L, 9007199254740992.0L);
+    TEST_f_f (rint, 9007199254740993.0000000000001L, 9007199254740993.0L);
+    TEST_f_f (rint, 9007199254740991.5000000000001L, 9007199254740991.0L);
+    TEST_f_f (rint, 9007199254740992.5000000000001L, 9007199254740992.0L);
+    TEST_f_f (rint, 9007199254740993.5000000000001L, 9007199254740993.0L);
+
+    TEST_f_f (rint, -9007199254740991.0000000000001L, -9007199254740992.0L);
+    TEST_f_f (rint, -9007199254740992.0000000000001L, -9007199254740993.0L);
+    TEST_f_f (rint, -9007199254740993.0000000000001L, -9007199254740994.0L);
+    TEST_f_f (rint, -9007199254740991.5000000000001L, -9007199254740992.0L);
+    TEST_f_f (rint, -9007199254740992.5000000000001L, -9007199254740993.0L);
+    TEST_f_f (rint, -9007199254740993.5000000000001L, -9007199254740994.0L);
+#endif
   }
 
   fesetround(save_round_mode);
@@ -4120,6 +4313,45 @@
     TEST_f_f (rint, -1.0, -1.0);
     TEST_f_f (rint, -1.5, -1.0);
     TEST_f_f (rint, -2.0, -2.0);
+    TEST_f_f (rint, 0.1, 1.0);
+    TEST_f_f (rint, 0.25, 1.0);
+    TEST_f_f (rint, 0.625, 1.0);
+    TEST_f_f (rint, -0.1, -0.0);
+    TEST_f_f (rint, -0.25, -0.0);
+    TEST_f_f (rint, -0.625, -0.0);
+#ifdef TEST_LDOUBLE
+  /* The result can only be represented in long double.  */
+    TEST_f_f (rint, 4503599627370495.5L, 4503599627370496.0L);
+    TEST_f_f (rint, 4503599627370496.25L, 4503599627370497.0L);
+    TEST_f_f (rint, 4503599627370496.5L, 4503599627370497.0L);
+    TEST_f_f (rint, 4503599627370496.75L, 4503599627370497.0L);
+    TEST_f_f (rint, 4503599627370497.5L, 4503599627370498.0L);
+    TEST_f_f (rint, 4503599627370494.5000000000001L, 4503599627370495.0L);
+    TEST_f_f (rint, 4503599627370495.5000000000001L, 4503599627370496.0L);
+    TEST_f_f (rint, 4503599627370496.5000000000001L, 4503599627370497.0L);
+    TEST_f_f (rint, -4503599627370495.5L, -4503599627370495.0L);
+    TEST_f_f (rint, -4503599627370496.25L, -4503599627370496.0L);
+    TEST_f_f (rint, -4503599627370496.5L, -4503599627370496.0L);
+    TEST_f_f (rint, -4503599627370496.75L, -4503599627370496.0L);
+    TEST_f_f (rint, -4503599627370497.5L, -4503599627370497.0L);
+    TEST_f_f (rint, -4503599627370494.5000000000001L, -4503599627370494.0L);
+    TEST_f_f (rint, -4503599627370495.5000000000001L, -4503599627370495.0L);
+    TEST_f_f (rint, -4503599627370496.5000000000001L, -4503599627370496.0L);
+
+    TEST_f_f (rint, 9007199254740991.0000000000001L, 9007199254740992.0L);
+    TEST_f_f (rint, 9007199254740992.0000000000001L, 9007199254740993.0L);
+    TEST_f_f (rint, 9007199254740993.0000000000001L, 9007199254740994.0L);
+    TEST_f_f (rint, 9007199254740991.5000000000001L, 9007199254740992.0L);
+    TEST_f_f (rint, 9007199254740992.5000000000001L, 9007199254740993.0L);
+    TEST_f_f (rint, 9007199254740993.5000000000001L, 9007199254740994.0L);
+
+    TEST_f_f (rint, -9007199254740991.0000000000001L, -9007199254740991.0L);
+    TEST_f_f (rint, -9007199254740992.0000000000001L, -9007199254740992.0L);
+    TEST_f_f (rint, -9007199254740993.0000000000001L, -9007199254740993.0L);
+    TEST_f_f (rint, -9007199254740991.5000000000001L, -9007199254740991.0L);
+    TEST_f_f (rint, -9007199254740992.5000000000001L, -9007199254740992.0L);
+    TEST_f_f (rint, -9007199254740993.5000000000001L, -9007199254740993.0L);
+#endif
   }
 
   fesetround(save_round_mode);
@@ -4142,6 +4374,12 @@
   TEST_f_f (round, -0.8L, -1.0);
   TEST_f_f (round, 1.5, 2.0);
   TEST_f_f (round, -1.5, -2.0);
+  TEST_f_f (round, 0.1, 0.0);
+  TEST_f_f (round, 0.25, 0.0);
+  TEST_f_f (round, 0.625, 1.0);
+  TEST_f_f (round, -0.1, -0.0);
+  TEST_f_f (round, -0.25, -0.0);
+  TEST_f_f (round, -0.625, -1.0);
   TEST_f_f (round, 2097152.5, 2097153);
   TEST_f_f (round, -2097152.5, -2097153);
 
@@ -4151,13 +4389,19 @@
   TEST_f_f (round, 4503599627370496.25L, 4503599627370496.0L);
   TEST_f_f (round, 4503599627370496.5L, 4503599627370497.0L); 
   TEST_f_f (round, 4503599627370496.75L, 4503599627370497.0L);
-  TEST_f_f (round, 4503599627370497.5L, 4503599627370498.0L);  
+  TEST_f_f (round, 4503599627370497.5L, 4503599627370498.0L);
+  TEST_f_f (round, 4503599627370494.5000000000001L, 4503599627370495.0L);
+  TEST_f_f (round, 4503599627370495.5000000000001L, 4503599627370496.0L);
+  TEST_f_f (round, 4503599627370496.5000000000001L, 4503599627370497.0L);
 
   TEST_f_f (round, -4503599627370495.5L, -4503599627370496.0L); 
   TEST_f_f (round, -4503599627370496.25L, -4503599627370496.0L); 
   TEST_f_f (round, -4503599627370496.5L, -4503599627370497.0L);
   TEST_f_f (round, -4503599627370496.75L, -4503599627370497.0L); 
   TEST_f_f (round, -4503599627370497.5L, -4503599627370498.0L);
+  TEST_f_f (round, -4503599627370494.5000000000001L, -4503599627370495.0L);
+  TEST_f_f (round, -4503599627370495.5000000000001L, -4503599627370496.0L);
+  TEST_f_f (round, -4503599627370496.5000000000001L, -4503599627370497.0L);
 
   TEST_f_f (round, 9007199254740991.5L, 9007199254740992.0L);
   TEST_f_f (round, 9007199254740992.25L, 9007199254740992.0L);
@@ -4171,6 +4415,20 @@
   TEST_f_f (round, -9007199254740992.75L, -9007199254740993.0L);
   TEST_f_f (round, -9007199254740993.5L, -9007199254740994.0L);
 
+  TEST_f_f (round, 9007199254740991.0000000000001L, 9007199254740991.0L);
+  TEST_f_f (round, 9007199254740992.0000000000001L, 9007199254740992.0L);
+  TEST_f_f (round, 9007199254740993.0000000000001L, 9007199254740993.0L);
+  TEST_f_f (round, 9007199254740991.5000000000001L, 9007199254740992.0L);
+  TEST_f_f (round, 9007199254740992.5000000000001L, 9007199254740993.0L);
+  TEST_f_f (round, 9007199254740993.5000000000001L, 9007199254740994.0L);
+
+  TEST_f_f (round, -9007199254740991.0000000000001L, -9007199254740991.0L);
+  TEST_f_f (round, -9007199254740992.0000000000001L, -9007199254740992.0L);
+  TEST_f_f (round, -9007199254740993.0000000000001L, -9007199254740993.0L);
+  TEST_f_f (round, -9007199254740991.5000000000001L, -9007199254740992.0L);
+  TEST_f_f (round, -9007199254740992.5000000000001L, -9007199254740993.0L);
+  TEST_f_f (round, -9007199254740993.5000000000001L, -9007199254740994.0L);
+
   TEST_f_f (round, 72057594037927935.5L, 72057594037927936.0L);
   TEST_f_f (round, 72057594037927936.25L, 72057594037927936.0L);
   TEST_f_f (round, 72057594037927936.5L, 72057594037927937.0L);
@@ -4541,7 +4799,11 @@
 
   TEST_f_f (trunc, 0, 0);
   TEST_f_f (trunc, minus_zero, minus_zero);
+  TEST_f_f (trunc, 0.1, 0);
+  TEST_f_f (trunc, 0.25, 0);
   TEST_f_f (trunc, 0.625, 0);
+  TEST_f_f (trunc, -0.1, minus_zero);
+  TEST_f_f (trunc, -0.25, minus_zero);
   TEST_f_f (trunc, -0.625, minus_zero);
   TEST_f_f (trunc, 1, 1);
   TEST_f_f (trunc, -1, -1);
@@ -4565,11 +4827,19 @@
   TEST_f_f (trunc, 4503599627370496.75L, 4503599627370496.0L);
   TEST_f_f (trunc, 4503599627370497.5L, 4503599627370497.0L);
 
+  TEST_f_f (trunc, 4503599627370494.5000000000001L, 4503599627370494.0L);
+  TEST_f_f (trunc, 4503599627370495.5000000000001L, 4503599627370495.0L);
+  TEST_f_f (trunc, 4503599627370496.5000000000001L, 4503599627370496.0L);
+  
   TEST_f_f (trunc, -4503599627370495.5L, -4503599627370495.0L);
   TEST_f_f (trunc, -4503599627370496.25L, -4503599627370496.0L);
   TEST_f_f (trunc, -4503599627370496.5L, -4503599627370496.0L);
   TEST_f_f (trunc, -4503599627370496.75L, -4503599627370496.0L);
   TEST_f_f (trunc, -4503599627370497.5L, -4503599627370497.0L);
+  
+  TEST_f_f (trunc, -4503599627370494.5000000000001L, -4503599627370494.0L);
+  TEST_f_f (trunc, -4503599627370495.5000000000001L, -4503599627370495.0L);
+  TEST_f_f (trunc, -4503599627370496.5000000000001L, -4503599627370496.0L);
 
   TEST_f_f (trunc, 9007199254740991.5L, 9007199254740991.0L);
   TEST_f_f (trunc, 9007199254740992.25L, 9007199254740992.0L);
@@ -4577,12 +4847,26 @@
   TEST_f_f (trunc, 9007199254740992.75L, 9007199254740992.0L);
   TEST_f_f (trunc, 9007199254740993.5L, 9007199254740993.0L);
 
+  TEST_f_f (trunc, 9007199254740991.0000000000001L, 9007199254740991.0L);
+  TEST_f_f (trunc, 9007199254740992.0000000000001L, 9007199254740992.0L);
+  TEST_f_f (trunc, 9007199254740993.0000000000001L, 9007199254740993.0L);
+  TEST_f_f (trunc, 9007199254740991.5000000000001L, 9007199254740991.0L);
+  TEST_f_f (trunc, 9007199254740992.5000000000001L, 9007199254740992.0L);
+  TEST_f_f (trunc, 9007199254740993.5000000000001L, 9007199254740993.0L);
+
   TEST_f_f (trunc, -9007199254740991.5L, -9007199254740991.0L);
   TEST_f_f (trunc, -9007199254740992.25L, -9007199254740992.0L);
   TEST_f_f (trunc, -9007199254740992.5L, -9007199254740992.0L);
   TEST_f_f (trunc, -9007199254740992.75L, -9007199254740992.0L);
   TEST_f_f (trunc, -9007199254740993.5L, -9007199254740993.0L);
 
+  TEST_f_f (trunc, -9007199254740991.0000000000001L, -9007199254740991.0L);
+  TEST_f_f (trunc, -9007199254740992.0000000000001L, -9007199254740992.0L);
+  TEST_f_f (trunc, -9007199254740993.0000000000001L, -9007199254740993.0L);
+  TEST_f_f (trunc, -9007199254740991.5000000000001L, -9007199254740991.0L);
+  TEST_f_f (trunc, -9007199254740992.5000000000001L, -9007199254740992.0L);
+  TEST_f_f (trunc, -9007199254740993.5000000000001L, -9007199254740993.0L);
+
   TEST_f_f (trunc, 72057594037927935.5L, 72057594037927935.0L);
   TEST_f_f (trunc, 72057594037927936.25L, 72057594037927936.0L);
   TEST_f_f (trunc, 72057594037927936.5L, 72057594037927936.0L);
diff -urN libc24-cvstip-20060301/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h libc24/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h
--- libc24-cvstip-20060301/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h	2006-01-27 18:07:25.000000000 -0600
+++ libc24/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h	2006-03-02 11:24:06.000000000 -0600
@@ -3,6 +3,7 @@
 #endif
 
 #include <sysdeps/ieee754/ldbl-128/math_ldbl.h>
+#include <ieee754.h>
 
 #define EXTRACT_IBM_EXTENDED_MANTISSA(hi64, lo64, expnt, ibm_ext_ldbl) \
   do									      \
@@ -122,3 +123,62 @@
       ibm_ext_ldbl = u.d;						      \
     }									      \
   while (0)
+  
+  
+/* Handy utility functions to pack/unpack/cononicalize and find the nearbyint
+   of long double implemented as double double.  */
+static inline long double
+__ldbl_pack (double a, double aa)
+{
+  union ibm_extended_long_double u;
+  u.dd[0] = a;
+  u.dd[1] = aa;
+  return u.d;
+}
+
+static inline void
+__ldbl_unpack (long double l, double *a, double *aa)
+{
+  union ibm_extended_long_double u;
+  u.d = l;
+  *a = u.dd[0];
+  *aa = u.dd[1];
+}
+
+
+/* Convert a finite long double to canonical form.
+   Does not handle +/-Inf properly.  */
+static inline void
+__ldbl_canonicalise (double *a, double *aa)
+{
+  double xh, xl;
+
+  xh = *a + *aa;
+  xl = (*a - xh) + *aa;
+  *a = xh;
+  *aa = xl;
+}
+
+/* Simple inline nearbyint (double) function .
+   Only works in the default rounding mode
+   but is useful in long double rounding functions.  */
+static inline double
+__ldbl_nearbyint (double a)
+{
+  double two52 = 0x10000000000000LL;
+
+  if (__builtin_expect ((__builtin_fabs (a) < two52), 1))
+    {
+      if (__builtin_expect ((a > 0.0), 1))
+	{
+	  a += two52;
+	  a -= two52;
+	}
+      else if (__builtin_expect ((a < 0.0), 1))
+	{
+	  a = two52 - a;
+	  a = -(a - two52);
+	}
+    }
+  return a;
+}
diff -urN libc24-cvstip-20060301/sysdeps/ieee754/ldbl-128ibm/s_ceill.c libc24/sysdeps/ieee754/ldbl-128ibm/s_ceill.c
--- libc24-cvstip-20060301/sysdeps/ieee754/ldbl-128ibm/s_ceill.c	2006-01-31 12:56:23.000000000 -0600
+++ libc24/sysdeps/ieee754/ldbl-128ibm/s_ceill.c	2006-03-01 17:05:50.000000000 -0600
@@ -19,7 +19,7 @@
    02111-1307 USA.  */
 
 #include <math.h>
-#include <fenv.h>
+#include <fenv_libc.h>
 #include <math_ldbl_opt.h>
 #include <float.h>
 #include <ieee754.h>
@@ -34,87 +34,58 @@
      long double x;
 #endif
 {
-  static const double TWO52 = 4503599627370496.0L;
-  int mode = fegetround();
-  union ibm_extended_long_double u;
+  double xh, xl, hi, lo;
 
-  u.d = x;
+  __ldbl_unpack (x, &xh, &xl);
 
-  if (fabs (u.dd[0]) < TWO52)
+  /* Return Inf, Nan, +/-0 unchanged.  */
+  if (__builtin_expect (xh != 0.0
+			&& __builtin_isless (__builtin_fabs (xh),
+					     __builtin_inf ()), 1))
     {
-      double high = u.dd[0];
-      fesetround(FE_UPWARD);
-      if (high > 0.0)
-	{
-	  high += TWO52;
-	  high -= TWO52;
-          if (high == -0.0) high = 0.0;
-	}
-      else if (high < 0.0)
-	{
-	  high -= TWO52;
-	  high += TWO52;
-          if (high == 0.0) high = -0.0;
-	}
-      u.dd[0] = high;
-      u.dd[1] = 0.0;
-      fesetround(mode);
-    }
-  else if (fabs (u.dd[1]) < TWO52 && u.dd[1] != 0.0)
-    {
-      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)
-	{
-	  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_UPWARD);
-	  low += TWO52;
-	  low -= TWO52;
-          fesetround(mode);
-	}
-      else if (u.dd[0] < 0.0)
-	{
-	  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_UPWARD);
-	  low -= TWO52;
-	  low += TWO52;
-          fesetround(mode);
-	}
-      u.dd[0] = high + low;
-      u.dd[1] = high - u.dd[0] + low;
+      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_canonicalise (&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 non-zero.  nearbyint
+         rounds values to the nearest integer, and values halfway
+         between integers to the nearest even integer.  ceill must
+         round towards +Inf.  */
+      xh -= lo;
+      __ldbl_canonicalise (&xh, &xl);
+
+      if (xh > 0.0 || (xh == 0.0 && xl > 0.0))
+	lo += 1.0;
+
+      /* 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_canonicalise (&xh, &xl);
+
+      /* Ensure we return -0 rather than +0 when appropriate.  */
+      if (orig_xh < 0.0)
+	xh = -__builtin_fabs (xh);
+
+      fesetround (save_round);
     }
 
-  return u.d;
+  return __ldbl_pack (xh, xl);
 }
 
 long_double_symbol (libm, __ceill, ceill);
diff -urN libc24-cvstip-20060301/sysdeps/ieee754/ldbl-128ibm/s_floorl.c libc24/sysdeps/ieee754/ldbl-128ibm/s_floorl.c
--- libc24-cvstip-20060301/sysdeps/ieee754/ldbl-128ibm/s_floorl.c	2006-01-31 12:56:23.000000000 -0600
+++ libc24/sysdeps/ieee754/ldbl-128ibm/s_floorl.c	2006-03-01 17:05:50.000000000 -0600
@@ -19,7 +19,7 @@
    02111-1307 USA.  */
 
 #include <math.h>
-#include <fenv.h>
+#include <fenv_libc.h>
 #include <math_ldbl_opt.h>
 #include <float.h>
 #include <ieee754.h>
@@ -34,86 +34,52 @@
      long double x;
 #endif
 {
-  static const double TWO52 = 4503599627370496.0L;
-  int mode = fegetround();
-  union ibm_extended_long_double u;
+  double xh, xl, hi, lo;
 
-  u.d = x;
+  __ldbl_unpack (x, &xh, &xl);
 
-  if (fabs (u.dd[0]) < TWO52)
+  /* Return Inf, Nan, +/-0 unchanged.  */
+  if (__builtin_expect (xh != 0.0
+			&& __builtin_isless (__builtin_fabs (xh),
+					     __builtin_inf ()), 1))
     {
-      double high = u.dd[0];
-      fesetround(FE_DOWNWARD);
-      if (high > 0.0)
-	{
-	  high += TWO52;
-	  high -= TWO52;
-          if (high == -0.0) high = 0.0;
-	}
-      else if (high < 0.0)
-	{
-	  high -= TWO52;
-	  high += TWO52;
-          if (high == 0.0) high = -0.0;
-	}
-      u.dd[0] = high;
-      u.dd[1] = 0.0;
-      fesetround(mode);
-    }
-  else if (fabs (u.dd[1]) < TWO52 && u.dd[1] != 0.0)
-    {
-      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)
-	{
-	  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_DOWNWARD);
-	  low += TWO52;
-	  low -= TWO52;
-	}
-      else if (u.dd[0] < 0.0)
-	{
-	  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_DOWNWARD);
-	  low -= TWO52;
-	  low += TWO52;
-	}
-      fesetround(mode);
-      u.dd[0] = high + low;
-      u.dd[1] = high - u.dd[0] + low;
+      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.  */
+      hi = __ldbl_nearbyint (xh);
+
+      /* Subtract integral high part from the value.  */
+      xh -= hi;
+      __ldbl_canonicalise (&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 non-zero.  nearbyint
+         rounds values to the nearest integer, and values halfway
+         between integers to the nearest even integer.  floorl must
+         round towards -Inf.  */
+      xh -= lo;
+      __ldbl_canonicalise (&xh, &xl);
+
+      if (xh < 0.0 || (xh == 0.0 && xl < 0.0))
+	lo += -1.0;
+
+      /* 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_canonicalise (&xh, &xl);
+
+      fesetround (save_round);
     }
 
-  return u.d;
+  return __ldbl_pack (xh, xl);
 }
 
 long_double_symbol (libm, __floorl, floorl);
diff -urN libc24-cvstip-20060301/sysdeps/ieee754/ldbl-128ibm/s_rintl.c libc24/sysdeps/ieee754/ldbl-128ibm/s_rintl.c
--- libc24-cvstip-20060301/sysdeps/ieee754/ldbl-128ibm/s_rintl.c	2006-01-27 18:07:25.000000000 -0600
+++ libc24/sysdeps/ieee754/ldbl-128ibm/s_rintl.c	2006-03-01 17:05:50.000000000 -0600
@@ -22,6 +22,7 @@
    when it's coded in C.  */
 
 #include <math.h>
+#include <fenv_libc.h>
 #include <math_ldbl_opt.h>
 #include <float.h>
 #include <ieee754.h>
@@ -36,84 +37,83 @@
      long double x;
 #endif
 {
-  static const long double TWO52 = 4503599627370496.0L;
-  union ibm_extended_long_double u;
-  u.d = x;
+  double xh, xl, hi, lo;
 
-  if (fabs (u.dd[0]) < TWO52)
-    {
-      double high = u.dd[0];
-      if (high > 0.0)
-	{
-	  high += TWO52;
-	  high -= TWO52;
-          if (high == -0.0) high = 0.0;
-	}
-      else if (high < 0.0)
-	{
-	  high -= TWO52;
-	  high += TWO52;
-          if (high == 0.0) high = -0.0;
-	}
-      u.dd[0] = high;
-      u.dd[1] = 0.0;
-    }
-  else if (fabs (u.dd[1]) < TWO52 && u.dd[1] != 0.0)
+  __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, tau;
-      /* 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)
-	{
-	  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.  */
-	      
-	      tau = nextafter (u.dd[0], 0.0);
-	      tau = (u.dd[0] - tau) * 2.0;
-	      high = u.dd[0] - tau;
-	      low = u.dd[1] + tau;
-	    }
-	  low += TWO52;
-	  low -= TWO52;
-	}
-      else 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.  If the low double
+	 happens to be exactly 0.5 or -0.5, you might think that this
+	 subtraction could result in an incorrect conversion.  For
+	 instance, subtracting an odd number would cause this function
+	 to round in the wrong direction.  However, if we have a
+	 canonical long double with the low double 0.5 or -0.5, then the
+	 high double must be even.  */
+      xh -= hi;
+      __ldbl_canonicalise (&xh, &xl);
+
+      /* Now convert the low double, adjusted for any remainder from the
+	 high double.  */
+      lo = __ldbl_nearbyint (xh);
+
+      xh -= lo;
+      __ldbl_canonicalise (&xh, &xl);
+
+      switch (save_round)
 	{
-	  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.  */
-	      tau = nextafter (u.dd[0], 0.0);
-	      tau = (u.dd[0] - tau) * 2.0;
-	      high = u.dd[0] - tau;
-	      low = u.dd[1] + tau;
-	    }
-	  low = TWO52 - low;
-	  low = -(low - TWO52);
+	case FE_TONEAREST:
+	  if (xl > 0.0 && xh == 0.5)
+	    lo += 1.0;
+	  else if (xl < 0.0 && -xh == 0.5)
+	    lo -= 1.0;
+	  break;
+
+	case FE_TOWARDZERO:
+	  if (orig_xh < 0.0)
+	    goto do_up;
+	  /* Fall thru */
+
+	case FE_DOWNWARD:
+	  if (xh < 0.0 || (xh == 0.0 && xl < 0.0))
+	    lo -= 1.0;
+	  break;
+
+	case FE_UPWARD:
+	do_up:
+	  if (xh > 0.0 || (xh == 0.0 && xl > 0.0))
+	    lo += 1.0;
+	  break;
 	}
-      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_canonicalise (&xh, &xl);
+
+      /* Ensure we return -0 rather than +0 when appropriate.  */
+      if (orig_xh < 0.0)
+	xh = -__builtin_fabs (xh);
+
+      fesetround (save_round);
     }
 
-  return u.d;
+  return __ldbl_pack (xh, xl);
 }
 
 long_double_symbol (libm, __rintl, rintl);
diff -urN libc24-cvstip-20060301/sysdeps/ieee754/ldbl-128ibm/s_roundl.c libc24/sysdeps/ieee754/ldbl-128ibm/s_roundl.c
--- libc24-cvstip-20060301/sysdeps/ieee754/ldbl-128ibm/s_roundl.c	2006-01-27 18:07:25.000000000 -0600
+++ libc24/sysdeps/ieee754/ldbl-128ibm/s_roundl.c	2006-03-01 17:05:50.000000000 -0600
@@ -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 @@
      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_canonicalise (&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_canonicalise (&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_canonicalise (&xh, &xl);
+
+      fesetround (save_round);
     }
-  return u.d;
+
+  return __ldbl_pack (xh, xl);
 }
 
 long_double_symbol (libm, __roundl, roundl);
diff -urN libc24-cvstip-20060301/sysdeps/ieee754/ldbl-128ibm/s_truncl.c libc24/sysdeps/ieee754/ldbl-128ibm/s_truncl.c
--- libc24-cvstip-20060301/sysdeps/ieee754/ldbl-128ibm/s_truncl.c	2006-01-27 18:07:25.000000000 -0600
+++ libc24/sysdeps/ieee754/ldbl-128ibm/s_truncl.c	2006-03-01 17:05:50.000000000 -0600
@@ -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,83 +37,66 @@
      long double x;
 #endif
 {
-  static const double TWO52 = 4503599627370496.0L;
-  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] += TWO52;
-	  u.dd[0] -= TWO52;
-	}
-      else if (u.dd[0] < 0.0)
-	{
-	  u.dd[0] = TWO52 - u.dd[0];
-	  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_canonicalise (&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 non-zero.  nearbyint
+         rounds values to the nearest integer, and values halfway
+         between integers to the nearest even integer.  floorl must
+         round towards -Inf.  */
+      xh -= lo;
+      __ldbl_canonicalise (&xh, &xl);
+
+      if (orig_xh < 0.0)
 	{
-	  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 += TWO52;
-	  low -= TWO52;
-          fesetround(mode);
+	  if (xh > 0.0 || (xh == 0.0 && xl > 0.0))
+	    lo += 1.0;
 	}
-      else if (u.dd[0] < 0.0)
+      else
 	{
-	  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 = TWO52 - low;
-	  low = -(low - TWO52);
-          fesetround(mode);
+	  if (xh < 0.0 || (xh == 0.0 && xl < 0.0))
+	    lo -= 1.0;
 	}
-      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_canonicalise (&xh, &xl);
+
+      /* Ensure we return -0 rather than +0 when appropriate.  */
+      if (orig_xh < 0.0)
+	xh = -__builtin_fabs (xh);
+
+      fesetround (save_round);
     }
 
-  return u.d;
+  return __ldbl_pack (xh, xl);
 }
 
 long_double_symbol (libm, __truncl, truncl);
diff -urN libc24-cvstip-20060301/sysdeps/powerpc/fpu/fenv_libc.h libc24/sysdeps/powerpc/fpu/fenv_libc.h
--- libc24-cvstip-20060301/sysdeps/powerpc/fpu/fenv_libc.h	2001-07-05 23:56:02.000000000 -0500
+++ libc24/sysdeps/powerpc/fpu/fenv_libc.h	2006-03-01 17:06:39.000000000 -0600
@@ -1,5 +1,5 @@
 /* Internal libc stuff for floating point environment routines.
-   Copyright (C) 1997 Free Software Foundation, Inc.
+   Copyright (C) 1997, 2006 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
    The GNU C Library is free software; you can redistribute it and/or
@@ -54,6 +54,42 @@
   unsigned int l[2];
 } fenv_union_t;
 
+
+static inline int
+___fegetround (void)
+{
+  int result;
+  asm ("mcrfs 7,7 ; mfcr %0" : "=r"(result) : : "cr7"); \
+  return result & 3;
+}
+#define fegetround ___fegetround
+
+static inline int
+___fesetround (int round)
+{
+  if ((unsigned int) round > 3)
+    return 1;
+  else if ((unsigned int) round < 2)
+    {
+       asm volatile ("mtfsb0 30" : : );
+       if ((unsigned int) round == 0)
+         asm volatile ("mtfsb0 31" : : );
+       else
+         asm volatile ("mtfsb1 31" : : );
+    }
+  else
+    {
+       asm volatile ("mtfsb1 30" : : );
+       if ((unsigned int) round == 2)
+         asm volatile ("mtfsb0 31" : : );
+       else
+         asm volatile ("mtfsb1 31" : : );
+    }
+
+  return 0;
+}
+#define fesetround ___fesetround
+
 /* Definitions of all the FPSCR bit numbers */
 enum {
   FPSCR_FX = 0,    /* exception summary */
diff -urN libc24-cvstip-20060301/sysdeps/powerpc/fpu/fesetround.c libc24/sysdeps/powerpc/fpu/fesetround.c
--- libc24-cvstip-20060301/sysdeps/powerpc/fpu/fesetround.c	2005-07-08 13:52:46.000000000 -0500
+++ libc24/sysdeps/powerpc/fpu/fesetround.c	2006-03-01 17:07:31.000000000 -0600
@@ -1,5 +1,5 @@
 /* Set current rounding direction.
-   Copyright (C) 1997, 2005 Free Software Foundation, Inc.
+   Copyright (C) 1997, 2005, 2006 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
@@ -20,22 +20,28 @@
 
 #include <fenv_libc.h>
 
+#undef fesetround
 int
 fesetround (int round)
 {
-  fenv_union_t u;
-
   if ((unsigned int) round > 3)
     return 1;
-
-  /* Get the current state.  */
-  u.fenv = fegetenv_register ();
-
-  /* Set the relevant bits.  */
-  u.l[1] = (u.l[1] & ~3)  |  (round & 3);
-
-  /* Put the new state in effect.  */
-  fesetenv_register (u.fenv);
+  else if ((unsigned int) round < 2)
+    {
+       asm volatile ("mtfsb0 30" : : );
+       if ((unsigned int) round == 0)
+         asm volatile ("mtfsb0 31" : : );
+       else
+         asm volatile ("mtfsb1 31" : : );
+    }
+  else
+    {
+       asm volatile ("mtfsb1 30" : : );
+       if ((unsigned int) round == 2)
+         asm volatile ("mtfsb0 31" : : );
+       else
+         asm volatile ("mtfsb1 31" : : );
+    }
 
   return 0;
 }
diff -urN libc24-cvstip-20060301/sysdeps/powerpc/fpu/math_ldbl.h libc24/sysdeps/powerpc/fpu/math_ldbl.h
--- libc24-cvstip-20060301/sysdeps/powerpc/fpu/math_ldbl.h	Wed Dec 31 18:00:00 1969
+++ libc24/sysdeps/powerpc/fpu/math_ldbl.h	Thu Mar 02 11:18:40 2006
@@ -0,0 +1,192 @@
+#ifndef _MATH_PRIVATE_H_
+#error "Never use <math_ldbl.h> directly; include <math_private.h> instead."
+#endif
+
+#include <sysdeps/ieee754/ldbl-128/math_ldbl.h>
+#include <ieee754.h>
+
+#define EXTRACT_IBM_EXTENDED_MANTISSA(hi64, lo64, expnt, ibm_ext_ldbl) \
+  do									      \
+    {									      \
+      /* We have 105 bits of mantissa plus one implicit digit.  Since	      \
+	 106 bits are representable without the rest using hexadecimal	      \
+	 digits we use only the implicit digits for the number before	      \
+	 the decimal point.  */						      \
+      unsigned long long hi, lo;					      \
+      int ediff;							      \
+      union ibm_extended_long_double eldbl;				      \
+      eldbl.d = ibm_ext_ldbl;						      \
+      expnt = eldbl.ieee.exponent - IBM_EXTENDED_LONG_DOUBLE_BIAS;	      \
+									      \
+      lo = ((long long)eldbl.ieee.mantissa2 << 32) | eldbl.ieee.mantissa3;    \
+      hi = ((long long)eldbl.ieee.mantissa0 << 32) | eldbl.ieee.mantissa1;    \
+      /* If the lower double is not a denomal or zero then set the hidden     \
+	 53rd bit.  */							      \
+      if (eldbl.ieee.exponent2 > 0x001)					      \
+	{								      \
+	  lo |= (1ULL << 52);						      \
+	  lo = lo << 7; /* pre-shift lo to match ieee854.  */		      \
+          /* The lower double is normalized separately from the upper.  We    \
+	     may need to adjust the lower manitissa to reflect this.  */      \
+	  ediff = eldbl.ieee.exponent - eldbl.ieee.exponent2;		      \
+	  if (ediff > 53)						      \
+	    lo = lo >> (ediff-53);					      \
+	}								      \
+      hi |= (1ULL << 52);						      \
+  									      \
+      if ((eldbl.ieee.negative != eldbl.ieee.negative2)			      \
+	  && ((eldbl.ieee.exponent2 != 0) && (lo != 0LL)))		      \
+	{								      \
+	  hi--;								      \
+	  lo = (1ULL << 60) - lo;					      \
+	  if (hi < (1ULL << 52))					      \
+	    {								      \
+	      /* we have a borrow from the hidden bit, so shift left 1.  */   \
+	      hi = (hi << 1) | (lo >> 59);				      \
+	      lo = 0xfffffffffffffffLL & (lo << 1);			      \
+	      expnt--;							      \
+	    }								      \
+	}								      \
+      lo64 = (hi << 60) | lo;						      \
+      hi64 = hi >> 4;							      \
+    }									      \
+  while (0)
+
+#define INSERT_IBM_EXTENDED_MANTISSA(ibm_ext_ldbl, sign, expnt, hi64, lo64) \
+  do									      \
+    {									      \
+      union ibm_extended_long_double u;					      \
+      unsigned long hidden2, lzcount;					      \
+      unsigned long long hi, lo;					      \
+									      \
+      u.ieee.negative = sign;						      \
+      u.ieee.negative2 = sign;						      \
+      u.ieee.exponent = expnt + IBM_EXTENDED_LONG_DOUBLE_BIAS;		      \
+      u.ieee.exponent2 = expnt-53 + IBM_EXTENDED_LONG_DOUBLE_BIAS;	      \
+      /* Expect 113 bits (112 bits + hidden) right justified in two longs.    \
+	 The low order 53 bits (52 + hidden) go into the lower double */      \
+      lo = (lo64 >> 7)& ((1ULL << 53) - 1);				      \
+      hidden2 = (lo64 >> 59) &  1ULL;					      \
+      /* The high order 53 bits (52 + hidden) go into the upper double */     \
+      hi = (lo64 >> 60) & ((1ULL << 11) - 1);				      \
+      hi |= (hi64 << 4);						      \
+  									      \
+      if (lo != 0LL)							      \
+	{								      \
+	  /* hidden2 bit of low double controls rounding of the high double.  \
+	     If hidden2 is '1' then round up hi and adjust lo (2nd mantissa)  \
+	     plus change the sign of the low double to compensate.  */	      \
+	  if (hidden2)							      \
+	    {								      \
+	      hi++;    							      \
+	      u.ieee.negative2 = !sign;					      \
+	      lo = (1ULL << 53) - lo;					      \
+	    }								      \
+	  /* The hidden bit of the lo mantissa is zero so we need to	      \
+	     normalize the it for the low double.  Shift it left until the    \
+	     hidden bit is '1' then adjust the 2nd exponent accordingly.  */  \
+									      \
+	  if (sizeof (lo) == sizeof (long))				      \
+	    lzcount = __builtin_clzl (lo);				      \
+	  else if ((lo >> 32) != 0)					      \
+	    lzcount = __builtin_clzl ((long) (lo >> 32));		      \
+	  else								      \
+	    lzcount = __builtin_clzl ((long) lo) + 32;			      \
+	  lzcount = lzcount - 11;					      \
+	  if (lzcount > 0)						      \
+	    {								      \
+	      int expnt2 = u.ieee.exponent2 - lzcount;			      \
+	      if (expnt2 >= 1)						      \
+		{							      \
+		  /* Not denormal.  Normalize and set low exponent.  */	      \
+		  lo = lo << lzcount; 					      \
+		  u.ieee.exponent2 = expnt2;				      \
+		}							      \
+	      else							      \
+		{							      \
+		  /* Is denormal.  */					      \
+		  lo = lo << (lzcount + expnt2);			      \
+		  u.ieee.exponent2 = 0;					      \
+		}							      \
+	    }								      \
+	}								      \
+      else 								      \
+	{								      \
+	  u.ieee.negative2 = 0;						      \
+	  u.ieee.exponent2 = 0;						      \
+	}								      \
+  									      \
+      u.ieee.mantissa3 = lo & ((1ULL << 32) - 1);			      \
+      u.ieee.mantissa2 = (lo >> 32) & ((1ULL << 20) - 1);		      \
+      u.ieee.mantissa1 = hi & ((1ULL << 32) - 1);			      \
+      u.ieee.mantissa0 = (hi >> 32) & ((1ULL << 20) - 1);		      \
+      ibm_ext_ldbl = u.d;						      \
+    }									      \
+  while (0)
+  
+  
+/* gcc generates disgusting code to pack and unpack long doubles.
+   This tells gcc that pack/unpack is really a nop.  We use fr1/fr2
+   because those are the regs used to pass/return a single
+   long double arg.  */
+static inline long double
+__ldbl_pack (double a, double aa)
+{
+  register long double x __asm__ ("fr1");
+  register double xh __asm__ ("fr1");
+  register double xl __asm__ ("fr2");
+  xh = a;
+  xl = aa;
+  __asm__ ("" : "=f" (x) : "f" (xh), "f" (xl));
+  return x;
+}
+
+static inline void
+__ldbl_unpack (long double l, double *a, double *aa)
+{
+  register long double x __asm__ ("fr1");
+  register double xh __asm__ ("fr1");
+  register double xl __asm__ ("fr2");
+  x = l;
+  __asm__ ("" : "=f" (xh), "=f" (xl) : "f" (x));
+  *a = xh;
+  *aa = xl;
+}
+
+
+/* Convert a finite long double to canonical form.
+   Does not handle +/-Inf properly.  */
+static inline void
+__ldbl_canonicalise (double *a, double *aa)
+{
+  double xh, xl;
+
+  xh = *a + *aa;
+  xl = (*a - xh) + *aa;
+  *a = xh;
+  *aa = xl;
+}
+
+/* Simple inline nearbyint (double) function .
+   Only works in the default rounding mode
+   but is useful in long double rounding functions.  */
+static inline double
+__ldbl_nearbyint (double a)
+{
+  double two52 = 0x10000000000000LL;
+
+  if (__builtin_expect ((__builtin_fabs (a) < two52), 1))
+    {
+      if (__builtin_expect ((a > 0.0), 1))
+	{
+	  a += two52;
+	  a -= two52;
+	}
+      else if (__builtin_expect ((a < 0.0), 1))
+	{
+	  a = two52 - a;
+	  a = -(a - two52);
+	}
+    }
+  return a;
+}
diff -urN libc24-cvstip-20060301/sysdeps/powerpc/powerpc64/fpu/s_rintl.S libc24/sysdeps/powerpc/powerpc64/fpu/s_rintl.S
--- libc24-cvstip-20060301/sysdeps/powerpc/powerpc64/fpu/s_rintl.S	Fri Jan 27 18:07:31 2006
+++ libc24/sysdeps/powerpc/powerpc64/fpu/s_rintl.S	Wed Dec 31 18:00:00 1969
@@ -1,113 +0,0 @@
-/* Round to int long double floating-point values.
-   IBM extended format long double version.
-   Copyright (C) 2004, 2006 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-
-   The GNU C Library is free software; you can redistribute it and/or
-   modify it under the terms of the GNU Lesser General Public
-   License as published by the Free Software Foundation; either
-   version 2.1 of the License, or (at your option) any later version.
-
-   The GNU C Library is distributed in the hope that it will be useful,
-   but WITHOUT ANY WARRANTY; without even the implied warranty of
-   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-   Lesser General Public License for more details.
-
-   You should have received a copy of the GNU Lesser General Public
-   License along with the GNU C Library; if not, write to the Free
-   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
-   02111-1307 USA.  */
-
-/* This has been coded in assembler because GCC makes such a mess of it
-   when it's coded in C.  */
-
-#include <sysdep.h>
-#include <math_ldbl_opt.h>
-
-	.section	".toc","aw"
-.LC0:	/* 2**52 */
-	.tc FD_43300000_0[TC],0x4330000000000000
-	.section	".text"
-
-ENTRY (__rintl)
-	lfd	fp13,.LC0@toc(2)
-	fabs	fp0,fp1
-	fsub	fp12,fp13,fp13	/* generate 0.0  */
-	fabs	fp9,fp2
-	fcmpu	cr7,fp0,fp13	/* if (fabs(x) > TWO52)  */
-	fcmpu	cr6,fp1,fp12	/* if (x > 0.0)  */
-	bnl-	cr7,.L2
-	fmr	fp2,fp12
-	bng-	cr6,.L1
-	fadd	fp1,fp1,fp13	/* x+= TWO52;  */
-	fsub	fp1,fp1,fp13	/* x-= TWO52;  */
-	fabs	fp1,fp1		/* if (x == 0.0)  */
-	blr			/* x = 0.0; */
-.L1:
-	bnllr-	cr6		/* if (x < 0.0)  */
-	fsub	fp1,fp1,fp13	/* x-= TWO52;  */
-	fadd	fp1,fp1,fp13	/* x+= TWO52;  */
-	fnabs	fp1,fp1		/* if (x == 0.0)  */
-	blr			/* x = -0.0; */
-
-/* The high double is > TWO52 so we need to round the low double and
-   perhaps the high double.  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.This gets a bit tricky so we use the
-   following algorithm:
-
-   tau = floor(x_high/TWO52);
-   x0 = x_high - tau;
-   x1 = x_low + tau;
-   r1 = rint(x1);
-   y_high = x0 + r1;
-   y_low = x0 - y_high + r1;
-   return y;  */
-.L2:
-	fcmpu	cr7,fp9,fp13	/* if (|x_low| > TWO52)  */
-	fcmpu	cr0,fp9,fp12	/* || (|x_low| == 0.0)  */
-	fcmpu	cr5,fp2,fp12	/* if (x_low > 0.0)  */
-	bgelr-	cr7		/*   return x;	*/
-	beqlr-  cr0
-	fdiv	fp8,fp1,fp13	/* x_high/TWO52  */
-	
-	bng-	cr6,.L6		/* if (x > 0.0)  */
-	fctidz	fp0,fp8
-	fcfid	fp8,fp0		/* tau = floor(x_high/TWO52);  */
-	fadd	fp8,fp8,fp8	/* tau++; Make tau even  */
-	bng	cr5,.L4		/* if (x_low > 0.0)  */
-	fmr	fp3,fp1
-	fmr	fp4,fp2
-	b	.L5
-.L4:				/* if (x_low < 0.0)  */
-	fsub	fp3,fp1,fp8	/* x0 = x_high - tau;  */
-	fadd	fp4,fp2,fp8	/* x1 = x_low + tau;  */
-.L5:
-	fadd	fp5,fp4,fp13	/* r1 = x1 + TWO52;  */
-	fsub	fp5,fp5,fp13	/* r1 = r1 - TWO52;  */
-	b	.L9
-.L6:				/* if (x < 0.0)  */
-	fctidz	fp0,fp8
-	fcfid	fp8,fp0		/* tau = floor(x_high/TWO52);  */
-	fadd	fp8,fp8,fp8	/* tau++; Make tau even  */	
-	bnl	cr5,.L7		/* if (x_low < 0.0)  */
-	fmr	fp3,fp1
-	fmr	fp4,fp2
-	b	.L8
-.L7:				/* if (x_low > 0.0)  */
-	fsub	fp3,fp1,fp8	/* x0 = x_high - tau;  */
-	fadd	fp4,fp2,fp8	/* x1 = x_low + tau;  */
-.L8:
-	fsub	fp5,fp13,fp4	/* r1 = TWO52 - x1;  */
-	fsub	fp0,fp5,fp13	/* r1 = - (r1 - TWO52);  */
-	fneg	fp5,fp0
-.L9:
-	fadd	fp1,fp3,fp5	/* y_high = x0 + r1;  */
-	fsub	fp2,fp3,fp1	/* y_low = x0 - y_high + r1;  */
-	fadd	fp2,fp2,fp5
-	blr
-END (__rintl)
-
-long_double_symbol (libm, __rintl, rintl)

Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]