lib/glibc
1
0
Fork 0
mirror of https://sourceware.org/git/glibc.git synced 2025-12-24 17:51:17 +03:00
Code Activity

[BZ #2423]

Browse Source 
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.
This commit is contained in:
Roland McGrath
2006-03-16 11:47:24 +00:00
parent 671ca699dd
commit 5c68d40169
15 changed files with 1046 additions and 622 deletions
Download Patch File Download Diff File Expand all files Collapse all files

39
ChangeLog
View File

@@ -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

334
math/libm-test.inc
View File

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

8
sysdeps/ieee754/ldbl-128ibm/e_fmodl.c
View File

@@ -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 */

5
sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c
View File

@@ -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

291
sysdeps/ieee754/ldbl-128ibm/math_ldbl.h
View File

@@ -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;
}

125
sysdeps/ieee754/ldbl-128ibm/s_ceill.c
View File

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

118
sysdeps/ieee754/ldbl-128ibm/s_floorl.c
View File

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

144
sysdeps/ieee754/ldbl-128ibm/s_rintl.c
View File

@@ -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)
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 = u.dd[0];
if (high > 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)
{
high += TWO52;
high -= TWO52;
if (high == -0.0) high = 0.0;
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;
}
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)
{
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)
{
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);
}
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);

122
sysdeps/ieee754/ldbl-128ibm/s_roundl.c
View File

@@ -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;
double xh, xl, hi, lo;
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)
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);

121
sysdeps/ieee754/ldbl-128ibm/s_truncl.c
View File

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

4
sysdeps/powerpc/fpu/fegetround.c
View File

@@ -23,7 +23,5 @@
int
fegetround (void)
{
int result;
asm ("mcrfs 7,7 ; mfcr %0" : "=r"(result) : : "cr7"); \
return result & 3;
return __fegetround();
}

37
sysdeps/powerpc/fpu/fenv_libc.h
View File

@@ -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 */

18
sysdeps/powerpc/fpu/fesetround.c
View File

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

189
sysdeps/powerpc/fpu/math_ldbl.h Normal file
View File

@@ -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;
}

113
sysdeps/powerpc/powerpc64/fpu/s_rintl.S
View File

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