2006-03-03 Steven Munroe <sjmunroe@us.ibm.com>

	    Alan Modra  <amodra@bigpond.net.au>

	[BZ #2423]
	* 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/fegetround.c: Use __fegetround.
	* sysdeps/powerpc/fpu/fesetround.c: Use __fesetround.

	* 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): Define inline utility
	functions for IBM long double format.
	* 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): Defined.

	* 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 --git a/ChangeLog b/ChangeLog
index 11e4518805b80a5b29ddc624a4b1406dd825d9a1..8e43a8cb852ff0dae3ea9fc6bd13d742dd00934e 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,37 @@
+2006-03-03  Steven Munroe  <sjmunroe@us.ibm.com>
+           Alan Modra  <amodra@bigpond.net.au>
+
+       [BZ #2423]
+       * 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/fegetround.c: Use __fegetround.
+       * sysdeps/powerpc/fpu/fesetround.c: Use __fesetround.
+
+       * 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): Define inline utility
+       functions for IBM long double format.
+       * 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): Defined.
+
+       * 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.
+
 2006-03-06  Roland McGrath  <roland@redhat.com>
 
        * version.h (VERSION): 2.4

diff --git a/math/libm-test.inc b/math/libm-test.inc
index 953459672549c92da45f02b343442bceacb69b72..e2a1d90bf651f8852fdbce5b9d9432f447cb2d7b 100644 (file)
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -1628,8 +1628,12 @@ ceil_test (void)
 
   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 @@ ceil_test (void)
   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 @@ ceil_test (void)
   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 @@ floor_test (void)
   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 @@ floor_test (void)
   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 @@ floor_test (void)
   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 @@ rint_test (void)
   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 @@ rint_test (void)
   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 @@ rint_test_tonearest (void)
     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 @@ rint_test_towardzero (void)
     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 @@ rint_test_downward (void)
     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 @@ rint_test_upward (void)
     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 @@ round_test (void)
   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 @@ round_test (void)
   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 @@ round_test (void)
   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 @@ trunc_test (void)
 
   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 @@ trunc_test (void)
   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 @@ trunc_test (void)
   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 --git a/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c b/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c
index 96668735b10d0261f4b79f2f9f3e035ebe4c84fa..e99b0bac34bcf271e88acae3fedcb2a3d24d6c2a 100644 (file)
--- a/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c
@@ -76,8 +76,8 @@ static long double one = 1.0, Zero[] = {0.0, -0.0,};
     /* Make the IBM extended format 105 bit mantissa look like the ieee854 112
        bit mantissa so the following operatations will give the correct
        result.  */
-        EXTRACT_IBM_EXTENDED_MANTISSA(hx, lx, temp, x);
-        EXTRACT_IBM_EXTENDED_MANTISSA(hy, ly, temp, y);
+        ldbl_extract_mantissa(&hx, &lx, &temp, x);
+        ldbl_extract_mantissa(&hy, &ly, &temp, y);
 
     /* set up {hx,lx}, {hy,ly} and align y to x */
        if(ix >= -1022)
@@ -127,7 +127,7 @@ static long double one = 1.0, Zero[] = {0.0, -0.0,};
            iy -= 1;
        }
        if(iy>= -1022) {        /* normalize output */
-           INSERT_IBM_EXTENDED_MANTISSA(x, (sx>>63), iy, hx, lx);
+           x = ldbl_insert_mantissa((sx>>63), iy, hx, lx);
        } else {                /* subnormal output */
            n = -1022 - iy;
            if(n<=48) {
@@ -138,7 +138,7 @@ static long double one = 1.0, Zero[] = {0.0, -0.0,};
            } else {
                lx = hx>>(n-64); hx = sx;
            }
-           INSERT_IBM_EXTENDED_MANTISSA(x, (sx>>63), iy, hx, lx);
+           x = ldbl_insert_mantissa((sx>>63), iy, hx, lx);
            x *= one;           /* create necessary signal */
        }
        return x;               /* exact output */

diff --git a/sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c b/sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c
index 416547c2208032f399494e9eda9106b640dd2e7c..8b1c976f3b038e7e4cd73a69ffbb314f58e24866 100644 (file)
--- a/sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c
@@ -199,7 +199,8 @@ int32_t __ieee754_rem_pio2l(long double x, long double *y)
 {
   long double z, w, t;
   double tx[8];
-  int64_t exp, n, ix, hx, ixd;
+  int exp;
+  int64_t n, ix, hx, ixd;
   u_int64_t lx, lxd;
 
   GET_LDOUBLE_WORDS64 (hx, lx, x);
@@ -243,7 +244,7 @@ int32_t __ieee754_rem_pio2l(long double x, long double *y)
      stored in a double array.  */
   /* Make the IBM extended format 105 bit mantissa look like the ieee854 112
      bit mantissa so the next operatation will give the correct result.  */
-  EXTRACT_IBM_EXTENDED_MANTISSA (ixd, lxd, exp, x);
+  ldbl_extract_mantissa (&ixd, &lxd, &exp, x);
   exp = exp - 23;
   /* This is faster than doing this in floating point, because we
      have to convert it to integers anyway and like this we can keep

diff --git a/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h b/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h
index e7e1b963e8ae7b66c2564d8c1ad6482ca4ce0d76..d055d6597ebe0cd480d82ff82203a9fc6cbd2d52 100644 (file)
--- a/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h
+++ b/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h
@@ -3,122 +3,179 @@
 #endif
 
 #include <sysdeps/ieee754/ldbl-128/math_ldbl.h>
+#include <ieee754.h>
+  
+static inline void
+ldbl_extract_mantissa (int64_t *hi64, u_int64_t *lo64, int *exp, long double x)
+{
+  /* We have 105 bits of mantissa plus one implicit digit.  Since
+     106 bits are representable we use the first implicit digit for
+     the number before the decimal point and the second implicit bit
+     as bit 53 of the mantissa.  */
+  unsigned long long hi, lo;
+  int ediff;
+  union ibm_extended_long_double eldbl;
+  eldbl.d = x;
+  *exp = eldbl.ieee.exponent - IBM_EXTENDED_LONG_DOUBLE_BIAS;
 
-#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)
+  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);
+         *exp = *exp - 1;
+       }
+    }
+  *lo64 = (hi << 60) | lo;
+  *hi64 = hi >> 4;
+}
 
-#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)
+static inline long double
+ldbl_insert_mantissa (int sign, int exp, int64_t hi64, u_int64_t lo64)
+{
+  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 = exp + IBM_EXTENDED_LONG_DOUBLE_BIAS;
+  u.ieee.exponent2 = exp-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);
+  return u.d;
+}
+  
+/* 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_canonicalize (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 --git a/sysdeps/ieee754/ldbl-128ibm/s_ceill.c b/sysdeps/ieee754/ldbl-128ibm/s_ceill.c
index a606548873e891a55ce110bab7571cb539e2a5ac..035e4f52ce07f82c676e35c66db20a17f4b3fd6d 100644 (file)
--- a/sysdeps/ieee754/ldbl-128ibm/s_ceill.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_ceill.c
@@ -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 @@ __ceill (x)
      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_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 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_canonicalize (&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_canonicalize (&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 --git a/sysdeps/ieee754/ldbl-128ibm/s_floorl.c b/sysdeps/ieee754/ldbl-128ibm/s_floorl.c
index 2be5e28698da9da547ba2fcd7ec6d76d7693afb4..4c4ae9b035fdada467a2c474715a1bf9ba4fa8a8 100644 (file)
--- a/sysdeps/ieee754/ldbl-128ibm/s_floorl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_floorl.c
@@ -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 @@ __floorl (x)
      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_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 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_canonicalize (&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_canonicalize (&xh, &xl);
+
+      fesetround (save_round);
     }
 
-  return u.d;
+  return ldbl_pack (xh, xl);
 }
 
 long_double_symbol (libm, __floorl, floorl);

diff --git a/sysdeps/ieee754/ldbl-128ibm/s_rintl.c b/sysdeps/ieee754/ldbl-128ibm/s_rintl.c
index 51b1fea2c40e12421667bf836e41ba720a2d2435..1f4e33f91cf3a0e9310970dc321c2fa181382d23 100644 (file)
--- a/sysdeps/ieee754/ldbl-128ibm/s_rintl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_rintl.c
@@ -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 @@ __rintl (x)
      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_canonicalize (&xh, &xl);
+
+      /* Now convert the low double, adjusted for any remainder from the
+        high double.  */
+      lo = ldbl_nearbyint (xh);
+
+      xh -= lo;
+      ldbl_canonicalize (&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_canonicalize (&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 --git a/sysdeps/ieee754/ldbl-128ibm/s_roundl.c b/sysdeps/ieee754/ldbl-128ibm/s_roundl.c
index 7fe84b87de1d6f01c3e7ea7c088de132c6ec4a18..0880e6ee227dcb86fba2642bce46850dcf51e198 100644 (file)
--- 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);

diff --git a/sysdeps/ieee754/ldbl-128ibm/s_truncl.c b/sysdeps/ieee754/ldbl-128ibm/s_truncl.c
index 14544af20661f924021944ffaf8925ba9768a040..d7bc47ee0de0cf5df75eb232ecbc2058f588ffa1 100644 (file)
--- a/sysdeps/ieee754/ldbl-128ibm/s_truncl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_truncl.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,83 +37,66 @@ __truncl (x)
      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)
-    {      
-      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)
+  /* 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 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_canonicalize (&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_canonicalize (&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 --git a/sysdeps/powerpc/fpu/fegetround.c b/sysdeps/powerpc/fpu/fegetround.c
index ae521fdf1090408604eb7abaaea0bcd6a726395e..93663ad5462b3dc3955e3d4ee763fa8e8cadc0ab 100644 (file)
--- a/sysdeps/powerpc/fpu/fegetround.c
+++ b/sysdeps/powerpc/fpu/fegetround.c
@@ -23,7 +23,5 @@
 int
 fegetround (void)
 {
-  int result;
-  asm ("mcrfs 7,7 ; mfcr %0" : "=r"(result) : : "cr7"); \
-  return result & 3;
+  return __fegetround();
 }

diff --git a/sysdeps/powerpc/fpu/fenv_libc.h b/sysdeps/powerpc/fpu/fenv_libc.h
index 7ae12a7d2b7b1f9bc57ad1ba69423cabb75e3190..fd5fc0c7675edcede25bdbb01dd4f8082697873d 100644 (file)
--- a/sysdeps/powerpc/fpu/fenv_libc.h
+++ b/sysdeps/powerpc/fpu/fenv_libc.h
@@ -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,41 @@ typedef union
   unsigned int l[2];
 } fenv_union_t;
 
+
+static inline int
+__fegetround (void)
+{
+  int result;
+  asm volatile ("mcrfs 7,7\n\t"
+               "mfcr  %0" : "=r"(result) : : "cr7");
+  return result & 3;
+}
+#define fegetround() __fegetround()
+
+static inline int
+__fesetround (int round)
+{
+  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(mode) __fesetround(mode)
+
 /* Definitions of all the FPSCR bit numbers */
 enum {
   FPSCR_FX = 0,    /* exception summary */

diff --git a/sysdeps/powerpc/fpu/fesetround.c b/sysdeps/powerpc/fpu/fesetround.c
index 67518d0df4f28952f1a8cc1cf03ab38f3c071f21..a7efa3bbb043bb1e797b5fa0267b0e876a0b241d 100644 (file)
--- a/sysdeps/powerpc/fpu/fesetround.c
+++ b/sysdeps/powerpc/fpu/fesetround.c
@@ -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,23 +20,13 @@
 
 #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);
-
-  return 0;
+  else
+    return __fesetround(round);
 }
 libm_hidden_def (fesetround)

diff --git a/sysdeps/powerpc/fpu/math_ldbl.h b/sysdeps/powerpc/fpu/math_ldbl.h
new file mode 100644 (file)
index 0000000..6cd6d0b
--- /dev/null
+++ b/sysdeps/powerpc/fpu/math_ldbl.h
@@ -0,0 +1,189 @@
+#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>
+  
+static inline void
+ldbl_extract_mantissa (int64_t *hi64, u_int64_t *lo64, int *exp, long double x)
+{
+  /* We have 105 bits of mantissa plus one implicit digit.  Since
+     106 bits are representable we use the first implicit digit for
+     the number before the decimal point and the second implicit bit
+     as bit 53 of the mantissa.  */
+  unsigned long long hi, lo;
+  int ediff;
+  union ibm_extended_long_double eldbl;
+  eldbl.d = x;
+  *exp = 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);
+         *exp = *exp - 1;
+       }
+    }
+  *lo64 = (hi << 60) | lo;
+  *hi64 = hi >> 4;
+}
+
+static inline long double
+ldbl_insert_mantissa (int sign, int exp, int64_t hi64, u_int64_t lo64)
+{
+  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 = exp + IBM_EXTENDED_LONG_DOUBLE_BIAS;
+  u.ieee.exponent2 = exp-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);
+  return u.d;
+}
+  
+/* 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_canonicalize (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;
+}