summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--ChangeLog39
-rw-r--r--math/libm-test.inc334
-rw-r--r--sysdeps/ieee754/ldbl-128ibm/e_fmodl.c8
-rw-r--r--sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c5
-rw-r--r--sysdeps/ieee754/ldbl-128ibm/math_ldbl.h291
-rw-r--r--sysdeps/ieee754/ldbl-128ibm/s_ceill.c125
-rw-r--r--sysdeps/ieee754/ldbl-128ibm/s_floorl.c118
-rw-r--r--sysdeps/ieee754/ldbl-128ibm/s_rintl.c144
-rw-r--r--sysdeps/ieee754/ldbl-128ibm/s_roundl.c124
-rw-r--r--sysdeps/ieee754/ldbl-128ibm/s_truncl.c121
-rw-r--r--sysdeps/powerpc/fpu/fegetround.c4
-rw-r--r--sysdeps/powerpc/fpu/fenv_libc.h37
-rw-r--r--sysdeps/powerpc/fpu/fesetround.c18
-rw-r--r--sysdeps/powerpc/fpu/math_ldbl.h189
-rw-r--r--sysdeps/powerpc/powerpc64/fpu/s_rintl.S113
15 files changed, 1047 insertions, 623 deletions
diff --git a/ChangeLog b/ChangeLog
index 5eb228aabf..a9ece8f8f7 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,42 @@
+2006-03-07 Jakub Jelinek <jakub@redhat.com>
+
+ [BZ #2423]
+ * math/libm-test.inc [TEST_LDOUBLE] (ceil_test, floor_test, rint_test,
+ round_test, trunc_test): Only run some of the new tests if
+ LDBL_MANT_DIG > 100.
+
+2006-03-03 Steven Munroe <sjmunroe@us.ibm.com>
+ Alan Modra <amodra@bigpond.net.au>
+
+ * sysdeps/powerpc/fpu/fenv_libc.h (__fegetround, __fesetround):
+ Define inline implementations.
+ * sysdeps/powerpc/fpu/fegetround.c: Use __fegetround.
+ * sysdeps/powerpc/fpu/fesetround.c: Use __fesetround.
+
+ * sysdeps/powerpc/fpu/math_ldbl.h: New file.
+
+ [BZ #2423]
+ * math/libm-test.inc [TEST_LDOUBLE] (ceil_test, floor_test, rint_test,
+ round_test, trunc_test): Add new tests.
+ * sysdeps/ieee754/ldbl-128ibm/math_ldbl.h
+ (EXTRACT_IBM_EXTENDED_MANTISSA, INSERT_IBM_EXTENDED_MANTISSA):
+ Removed, replaced with ...
+ (ldbl_extract_mantissa, ldbl_insert_mantissa, ldbl_pack, ldbl_unpack,
+ ldbl_canonicalise, ldbl_nearbyint): New functions.
+ * sysdeps/ieee754/ldbl-128ibm/e_fmodl.c (__ieee754_fmodl): Replace
+ EXTRACT_IBM_EXTENDED_MANTISSA and INSERT_IBM_EXTENDED_MANTISSA
+ with ldbl_extract_mantissa and ldbl_insert_mantissa.
+ * sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c (__ieee754_rem_pio2l):
+ Replace EXTRACT_IBM_EXTENDED_MANTISSA with ldbl_extract_mantissa.
+ (ldbl_extract_mantissa, ldbl_insert_mantissa): New inline functions.
+ * sysdeps/ieee754/ldbl-128ibm/s_ceill.c (__ceill): Handle rounding
+ that spans doubles in IBM long double format.
+ * sysdeps/ieee754/ldbl-128ibm/s_floorl.c: Likewise.
+ * sysdeps/ieee754/ldbl-128ibm/s_rintl.c: Likewise.
+ * sysdeps/ieee754/ldbl-128ibm/s_roundl.c: Likewise.
+ * sysdeps/ieee754/ldbl-128ibm/s_truncl.c: Likewise.
+ * sysdeps/powerpc/powerpc64/fpu/s_rintl.S: File removed.
+
2004-12-09 Randolph Chung <tausq@debian.org>
* sysdeps/unix/sysv/linux/kernel-features.h (__ASSUME_UTIMES): Don't
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 9534596725..b144796f2a 100644
--- 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. */
@@ -1645,6 +1649,15 @@ ceil_test (void)
TEST_f_f (ceil, -4503599627370496.75L, -4503599627370496.0L);
TEST_f_f (ceil, -4503599627370497.5L, -4503599627370497.0L);
+# if LDBL_MANT_DIG > 100
+ 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);
+# endif
+
TEST_f_f (ceil, 9007199254740991.5L, 9007199254740992.0L);
TEST_f_f (ceil, 9007199254740992.25L, 9007199254740993.0L);
TEST_f_f (ceil, 9007199254740992.5L, 9007199254740993.0L);
@@ -1657,6 +1670,22 @@ ceil_test (void)
TEST_f_f (ceil, -9007199254740992.75L, -9007199254740992.0L);
TEST_f_f (ceil, -9007199254740993.5L, -9007199254740993.0L);
+# if LDBL_MANT_DIG > 100
+ 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);
+# endif
+
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 +2657,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 +2671,22 @@ 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);
+# if LDBL_MANT_DIG > 100
+ TEST_f_f (floor, 4503599627370494.5000000000001L, 4503599627370494.0L);
+ TEST_f_f (floor, 4503599627370495.5000000000001L, 4503599627370495.0L);
+ TEST_f_f (floor, 4503599627370496.5000000000001L, 4503599627370496.0L);
+# endif
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);
+# if LDBL_MANT_DIG > 100
+ TEST_f_f (floor, -4503599627370494.5000000000001L, -4503599627370495.0L);
+ TEST_f_f (floor, -4503599627370495.5000000000001L, -4503599627370496.0L);
+ TEST_f_f (floor, -4503599627370496.5000000000001L, -4503599627370497.0L);
+# endif
TEST_f_f (floor, 9007199254740991.5L, 9007199254740991.0L);
TEST_f_f (floor, 9007199254740992.25L, 9007199254740992.0L);
@@ -2652,12 +2694,30 @@ floor_test (void)
TEST_f_f (floor, 9007199254740992.75L, 9007199254740992.0L);
TEST_f_f (floor, 9007199254740993.5L, 9007199254740993.0L);
+# if LDBL_MANT_DIG > 100
+ 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);
+# endif
+
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);
+# if LDBL_MANT_DIG > 100
+ 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);
+# endif
+
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 +4031,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 +4044,38 @@ 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);
+
+# if LDBL_MANT_DIG > 100
+ TEST_f_f (rint, 4503599627370494.5000000000001L, 4503599627370495.0L);
+ TEST_f_f (rint, 4503599627370495.5000000000001L, 4503599627370496.0L);
+ TEST_f_f (rint, 4503599627370496.5000000000001L, 4503599627370497.0L);
+# endif
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);
+
+# if LDBL_MANT_DIG > 100
+ 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
TEST_f_f (rint, 9007199254740991.5L, 9007199254740992.0L);
TEST_f_f (rint, 9007199254740992.25L, 9007199254740992.0L);
@@ -4039,6 +4131,49 @@ 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);
+# if LDBL_MANT_DIG > 100
+ TEST_f_f (rint, 4503599627370494.5000000000001L, 4503599627370495.0L);
+ TEST_f_f (rint, 4503599627370495.5000000000001L, 4503599627370496.0L);
+ TEST_f_f (rint, 4503599627370496.5000000000001L, 4503599627370497.0L);
+# endif
+ 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);
+# if LDBL_MANT_DIG > 100
+ 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
+#endif
}
fesetround(save_round_mode);
@@ -4066,6 +4201,49 @@ 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);
+# if LDBL_MANT_DIG > 100
+ TEST_f_f (rint, 4503599627370494.5000000000001L, 4503599627370494.0L);
+ TEST_f_f (rint, 4503599627370495.5000000000001L, 4503599627370495.0L);
+ TEST_f_f (rint, 4503599627370496.5000000000001L, 4503599627370496.0L);
+# endif
+ 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);
+# if LDBL_MANT_DIG > 100
+ 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
+#endif
}
fesetround(save_round_mode);
@@ -4093,6 +4271,49 @@ 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);
+# if LDBL_MANT_DIG > 100
+ TEST_f_f (rint, 4503599627370494.5000000000001L, 4503599627370494.0L);
+ TEST_f_f (rint, 4503599627370495.5000000000001L, 4503599627370495.0L);
+ TEST_f_f (rint, 4503599627370496.5000000000001L, 4503599627370496.0L);
+# endif
+ 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);
+# if LDBL_MANT_DIG > 100
+ 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
+#endif
}
fesetround(save_round_mode);
@@ -4120,6 +4341,49 @@ 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);
+# if LDBL_MANT_DIG > 100
+ TEST_f_f (rint, 4503599627370494.5000000000001L, 4503599627370495.0L);
+ TEST_f_f (rint, 4503599627370495.5000000000001L, 4503599627370496.0L);
+ TEST_f_f (rint, 4503599627370496.5000000000001L, 4503599627370497.0L);
+# endif
+ 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);
+# if LDBL_MANT_DIG > 100
+ 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
+#endif
}
fesetround(save_round_mode);
@@ -4142,6 +4406,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 +4421,23 @@ 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);
+# if LDBL_MANT_DIG > 100
+ TEST_f_f (round, 4503599627370494.5000000000001L, 4503599627370495.0L);
+ TEST_f_f (round, 4503599627370495.5000000000001L, 4503599627370496.0L);
+ TEST_f_f (round, 4503599627370496.5000000000001L, 4503599627370497.0L);
+# endif
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);
+# if LDBL_MANT_DIG > 100
+ TEST_f_f (round, -4503599627370494.5000000000001L, -4503599627370495.0L);
+ TEST_f_f (round, -4503599627370495.5000000000001L, -4503599627370496.0L);
+ TEST_f_f (round, -4503599627370496.5000000000001L, -4503599627370497.0L);
+# endif
TEST_f_f (round, 9007199254740991.5L, 9007199254740992.0L);
TEST_f_f (round, 9007199254740992.25L, 9007199254740992.0L);
@@ -4171,6 +4451,22 @@ round_test (void)
TEST_f_f (round, -9007199254740992.75L, -9007199254740993.0L);
TEST_f_f (round, -9007199254740993.5L, -9007199254740994.0L);
+# if LDBL_MANT_DIG > 100
+ 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);
+# endif
+
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 +4837,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 +4865,23 @@ trunc_test (void)
TEST_f_f (trunc, 4503599627370496.75L, 4503599627370496.0L);
TEST_f_f (trunc, 4503599627370497.5L, 4503599627370497.0L);
+# if LDBL_MANT_DIG > 100
+ TEST_f_f (trunc, 4503599627370494.5000000000001L, 4503599627370494.0L);
+ TEST_f_f (trunc, 4503599627370495.5000000000001L, 4503599627370495.0L);
+ TEST_f_f (trunc, 4503599627370496.5000000000001L, 4503599627370496.0L);
+# endif
+
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);
+
+# if LDBL_MANT_DIG > 100
+ TEST_f_f (trunc, -4503599627370494.5000000000001L, -4503599627370494.0L);
+ TEST_f_f (trunc, -4503599627370495.5000000000001L, -4503599627370495.0L);
+ TEST_f_f (trunc, -4503599627370496.5000000000001L, -4503599627370496.0L);
+# endif
TEST_f_f (trunc, 9007199254740991.5L, 9007199254740991.0L);
TEST_f_f (trunc, 9007199254740992.25L, 9007199254740992.0L);
@@ -4577,12 +4889,30 @@ trunc_test (void)
TEST_f_f (trunc, 9007199254740992.75L, 9007199254740992.0L);
TEST_f_f (trunc, 9007199254740993.5L, 9007199254740993.0L);
+# if LDBL_MANT_DIG > 100
+ 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);
+# endif
+
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);
+# if LDBL_MANT_DIG > 100
+ 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);
+# endif
+
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 96668735b1..e99b0bac34 100644
--- 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 416547c220..8b1c976f3b 100644
--- 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 e7e1b963e8..d055d6597e 100644
--- 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 a606548873..035e4f52ce 100644
--- 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 2be5e28698..4c4ae9b035 100644
--- 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 51b1fea2c4..1f4e33f91c 100644
--- 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 7fe84b87de..0880e6ee22 100644
--- a/sysdeps/ieee754/ldbl-128ibm/s_roundl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_roundl.c
@@ -22,7 +22,7 @@
when it's coded in C. */
#include <math.h>
-#include <fenv.h>
+#include <fenv_libc.h>
#include <math_ldbl_opt.h>
#include <float.h>
#include <ieee754.h>
@@ -37,84 +37,62 @@ __roundl (x)
long double x;
#endif
{
- static const double TWO52 = 4503599627370496.0;
- static const double HALF = 0.5;
- int mode = fegetround();
- union ibm_extended_long_double u;
- u.d = x;
-
- if (fabs (u.dd[0]) < TWO52)
- {
- fesetround(FE_TOWARDZERO);
- if (u.dd[0] > 0.0)
- {
- u.dd[0] += HALF;
- u.dd[0] += TWO52;
- u.dd[0] -= TWO52;
- }
- else if (u.dd[0] < 0.0)
- {
- u.dd[0] = TWO52 - (u.dd[0] - HALF);
- u.dd[0] = -(u.dd[0] - TWO52);
- }
- u.dd[1] = 0.0;
- fesetround(mode);
- }
- else if (fabs (u.dd[1]) < TWO52 && u.dd[1] != 0.0)
+ double xh, xl, hi, lo;
+
+ ldbl_unpack (x, &xh, &xl);
+
+ /* Return Inf, Nan, +/-0 unchanged. */
+ if (__builtin_expect (xh != 0.0
+ && __builtin_isless (__builtin_fabs (xh),
+ __builtin_inf ()), 1))
{
- double high, low;
- /* In this case we have to round the low double and handle any
- adjustment to the high double that may be caused by rounding
- (up). This is complicated by the fact that the high double
- may already be rounded and the low double may have the
- opposite sign to compensate. */
- if (u.dd[0] > 0.0)
+ double orig_xh;
+ int save_round = fegetround ();
+
+ /* Long double arithmetic, including the canonicalisation below,
+ only works in round-to-nearest mode. */
+ fesetround (FE_TONEAREST);
+
+ /* Convert the high double to integer. */
+ orig_xh = xh;
+ hi = ldbl_nearbyint (xh);
+
+ /* Subtract integral high part from the value. */
+ xh -= hi;
+ ldbl_canonicalize (&xh, &xl);
+
+ /* Now convert the low double, adjusted for any remainder from the
+ high double. */
+ lo = ldbl_nearbyint (xh);
+
+ /* Adjust the result when the remainder is exactly 0.5. nearbyint
+ rounds values halfway between integers to the nearest even
+ integer. roundl must round away from zero.
+ Also correct cases where nearbyint returns an incorrect value
+ for LO. */
+ xh -= lo;
+ ldbl_canonicalize (&xh, &xl);
+ if (xh == 0.5)
{
- if (u.dd[1] > 0.0)
- {
- /* If the high/low doubles are the same sign then simply
- round the low double. */
- high = u.dd[0];
- low = u.dd[1];
- }
- else if (u.dd[1] < 0.0)
- {
- /* Else the high double is pre rounded and we need to
- adjust for that. */
- high = nextafter (u.dd[0], 0.0);
- low = u.dd[1] + (u.dd[0] - high);
- }
- fesetround(FE_TOWARDZERO);
- low += HALF;
- low += TWO52;
- low -= TWO52;
+ if (xl > 0.0 || (xl == 0.0 && orig_xh > 0.0))
+ lo += 1.0;
}
- else if (u.dd[0] < 0.0)
+ else if (-xh == 0.5)
{
- if (u.dd[1] < 0.0)
- {
- /* If the high/low doubles are the same sign then simply
- round the low double. */
- high = u.dd[0];
- low = u.dd[1];
- }
- else if (u.dd[1] > 0.0)
- {
- /* Else the high double is pre rounded and we need to
- adjust for that. */
- high = nextafter (u.dd[0], 0.0);
- low = u.dd[1] + (u.dd[0] - high);
- }
- fesetround(FE_TOWARDZERO);
- low -= HALF;
- low = TWO52 - low;
- low = -(low - TWO52);
+ if (xl < 0.0 || (xl == 0.0 && orig_xh < 0.0))
+ lo -= 1.0;
}
- fesetround(mode);
- u.dd[0] = high + low;
- u.dd[1] = high - u.dd[0] + low;
+
+ /* Ensure the final value is canonical. In certain cases,
+ rounding causes hi,lo calculated so far to be non-canonical. */
+ xh = hi;
+ xl = lo;
+ ldbl_canonicalize (&xh, &xl);
+
+ fesetround (save_round);
}
- return u.d;
+
+ return ldbl_pack (xh, xl);
}
long_double_symbol (libm, __roundl, roundl);
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_truncl.c b/sysdeps/ieee754/ldbl-128ibm/s_truncl.c
index 14544af206..d7bc47ee0d 100644
--- 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 ae521fdf10..93663ad546 100644
--- 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 7ae12a7d2b..fd5fc0c767 100644
--- 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 67518d0df4..a7efa3bbb0 100644
--- 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
index 0000000000..6cd6d0bdfe
--- /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;
+}
diff --git a/sysdeps/powerpc/powerpc64/fpu/s_rintl.S b/sysdeps/powerpc/powerpc64/fpu/s_rintl.S
deleted file mode 100644
index 2ca2d4481f..0000000000
--- a/sysdeps/powerpc/powerpc64/fpu/s_rintl.S
+++ /dev/null
@@ -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)
View Lorry Controller Status