summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorJakub Jelinek <jakub@redhat.com>2006-03-07 08:24:12 +0000
committerJakub Jelinek <jakub@redhat.com>2006-03-07 08:24:12 +0000
commit4133f08e15452432e4ec0d384f8a7be1199901bb (patch)
tree4c85cfcf94198e172c0466a5c251b6a7399cdc56
parentd1e6d2178d0b242c9ad40f77216f2f0cae57b97a (diff)
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.
-rw-r--r--ChangeLog34
-rw-r--r--math/libm-test.inc288
-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
14 files changed, 996 insertions, 510 deletions
diff --git a/ChangeLog b/ChangeLog
index 11e4518805..8e43a8cb85 100644
--- 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 9534596725..e2a1d90bf6 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. */
@@ -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 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;
+}