diff options
| -rw-r--r-- | ChangeLog | 39 | ||||
| -rw-r--r-- | math/libm-test.inc | 334 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/e_fmodl.c | 8 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c | 5 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/math_ldbl.h | 291 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/s_ceill.c | 125 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/s_floorl.c | 118 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/s_rintl.c | 144 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/s_roundl.c | 124 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/s_truncl.c | 121 | ||||
| -rw-r--r-- | sysdeps/powerpc/fpu/fegetround.c | 4 | ||||
| -rw-r--r-- | sysdeps/powerpc/fpu/fenv_libc.h | 37 | ||||
| -rw-r--r-- | sysdeps/powerpc/fpu/fesetround.c | 18 | ||||
| -rw-r--r-- | sysdeps/powerpc/fpu/math_ldbl.h | 189 | ||||
| -rw-r--r-- | sysdeps/powerpc/powerpc64/fpu/s_rintl.S | 113 |
15 files changed, 1047 insertions, 623 deletions
@@ -1,3 +1,42 @@ +2006-03-07 Jakub Jelinek <jakub@redhat.com> + + [BZ #2423] + * math/libm-test.inc [TEST_LDOUBLE] (ceil_test, floor_test, rint_test, + round_test, trunc_test): Only run some of the new tests if + LDBL_MANT_DIG > 100. + +2006-03-03 Steven Munroe <sjmunroe@us.ibm.com> + Alan Modra <amodra@bigpond.net.au> + + * sysdeps/powerpc/fpu/fenv_libc.h (__fegetround, __fesetround): + Define inline implementations. + * sysdeps/powerpc/fpu/fegetround.c: Use __fegetround. + * sysdeps/powerpc/fpu/fesetround.c: Use __fesetround. + + * sysdeps/powerpc/fpu/math_ldbl.h: New file. + + [BZ #2423] + * math/libm-test.inc [TEST_LDOUBLE] (ceil_test, floor_test, rint_test, + round_test, trunc_test): Add new tests. + * sysdeps/ieee754/ldbl-128ibm/math_ldbl.h + (EXTRACT_IBM_EXTENDED_MANTISSA, INSERT_IBM_EXTENDED_MANTISSA): + Removed, replaced with ... + (ldbl_extract_mantissa, ldbl_insert_mantissa, ldbl_pack, ldbl_unpack, + ldbl_canonicalise, ldbl_nearbyint): New functions. + * sysdeps/ieee754/ldbl-128ibm/e_fmodl.c (__ieee754_fmodl): Replace + EXTRACT_IBM_EXTENDED_MANTISSA and INSERT_IBM_EXTENDED_MANTISSA + with ldbl_extract_mantissa and ldbl_insert_mantissa. + * sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c (__ieee754_rem_pio2l): + Replace EXTRACT_IBM_EXTENDED_MANTISSA with ldbl_extract_mantissa. + (ldbl_extract_mantissa, ldbl_insert_mantissa): New inline functions. + * sysdeps/ieee754/ldbl-128ibm/s_ceill.c (__ceill): Handle rounding + that spans doubles in IBM long double format. + * sysdeps/ieee754/ldbl-128ibm/s_floorl.c: Likewise. + * sysdeps/ieee754/ldbl-128ibm/s_rintl.c: Likewise. + * sysdeps/ieee754/ldbl-128ibm/s_roundl.c: Likewise. + * sysdeps/ieee754/ldbl-128ibm/s_truncl.c: Likewise. + * sysdeps/powerpc/powerpc64/fpu/s_rintl.S: File removed. + 2004-12-09 Randolph Chung <tausq@debian.org> * sysdeps/unix/sysv/linux/kernel-features.h (__ASSUME_UTIMES): Don't diff --git a/math/libm-test.inc b/math/libm-test.inc index 9534596..b144796 100644 --- a/math/libm-test.inc +++ b/math/libm-test.inc @@ -1628,8 +1628,12 @@ ceil_test (void) TEST_f_f (ceil, M_PIl, 4.0); TEST_f_f (ceil, -M_PIl, -3.0); + TEST_f_f (ceil, 0.1, 1.0); TEST_f_f (ceil, 0.25, 1.0); + TEST_f_f (ceil, 0.625, 1.0); + TEST_f_f (ceil, -0.1, minus_zero); TEST_f_f (ceil, -0.25, minus_zero); + TEST_f_f (ceil, -0.625, minus_zero); #ifdef TEST_LDOUBLE /* The result can only be represented in long double. */ @@ -1645,6 +1649,15 @@ ceil_test (void) TEST_f_f (ceil, -4503599627370496.75L, -4503599627370496.0L); TEST_f_f (ceil, -4503599627370497.5L, -4503599627370497.0L); +# if LDBL_MANT_DIG > 100 + TEST_f_f (ceil, 4503599627370494.5000000000001L, 4503599627370495.0L); + TEST_f_f (ceil, 4503599627370495.5000000000001L, 4503599627370496.0L); + TEST_f_f (ceil, 4503599627370496.5000000000001L, 4503599627370497.0L); + TEST_f_f (ceil, -4503599627370494.5000000000001L, -4503599627370494.0L); + TEST_f_f (ceil, -4503599627370495.5000000000001L, -4503599627370495.0L); + TEST_f_f (ceil, -4503599627370496.5000000000001L, -4503599627370496.0L); +# endif + TEST_f_f (ceil, 9007199254740991.5L, 9007199254740992.0L); TEST_f_f (ceil, 9007199254740992.25L, 9007199254740993.0L); TEST_f_f (ceil, 9007199254740992.5L, 9007199254740993.0L); @@ -1657,6 +1670,22 @@ ceil_test (void) TEST_f_f (ceil, -9007199254740992.75L, -9007199254740992.0L); TEST_f_f (ceil, -9007199254740993.5L, -9007199254740993.0L); +# if LDBL_MANT_DIG > 100 + TEST_f_f (ceil, 9007199254740991.0000000000001L, 9007199254740992.0L); + TEST_f_f (ceil, 9007199254740992.0000000000001L, 9007199254740993.0L); + TEST_f_f (ceil, 9007199254740993.0000000000001L, 9007199254740994.0L); + TEST_f_f (ceil, 9007199254740991.5000000000001L, 9007199254740992.0L); + TEST_f_f (ceil, 9007199254740992.5000000000001L, 9007199254740993.0L); + TEST_f_f (ceil, 9007199254740993.5000000000001L, 9007199254740994.0L); + + TEST_f_f (ceil, -9007199254740991.0000000000001L, -9007199254740991.0L); + TEST_f_f (ceil, -9007199254740992.0000000000001L, -9007199254740992.0L); + TEST_f_f (ceil, -9007199254740993.0000000000001L, -9007199254740993.0L); + TEST_f_f (ceil, -9007199254740991.5000000000001L, -9007199254740991.0L); + TEST_f_f (ceil, -9007199254740992.5000000000001L, -9007199254740992.0L); + TEST_f_f (ceil, -9007199254740993.5000000000001L, -9007199254740993.0L); +# endif + TEST_f_f (ceil, 72057594037927935.5L, 72057594037927936.0L); TEST_f_f (ceil, 72057594037927936.25L, 72057594037927937.0L); TEST_f_f (ceil, 72057594037927936.5L, 72057594037927937.0L); @@ -2628,9 +2657,12 @@ floor_test (void) TEST_f_f (floor, M_PIl, 3.0); TEST_f_f (floor, -M_PIl, -4.0); + TEST_f_f (floor, 0.1, 0.0); TEST_f_f (floor, 0.25, 0.0); + TEST_f_f (floor, 0.625, 0.0); + TEST_f_f (floor, -0.1, -1.0); TEST_f_f (floor, -0.25, -1.0); - + TEST_f_f (floor, -0.625, -1.0); #ifdef TEST_LDOUBLE /* The result can only be represented in long double. */ @@ -2639,12 +2671,22 @@ floor_test (void) TEST_f_f (floor, 4503599627370496.5L, 4503599627370496.0L); TEST_f_f (floor, 4503599627370496.75L, 4503599627370496.0L); TEST_f_f (floor, 4503599627370497.5L, 4503599627370497.0L); +# if LDBL_MANT_DIG > 100 + TEST_f_f (floor, 4503599627370494.5000000000001L, 4503599627370494.0L); + TEST_f_f (floor, 4503599627370495.5000000000001L, 4503599627370495.0L); + TEST_f_f (floor, 4503599627370496.5000000000001L, 4503599627370496.0L); +# endif TEST_f_f (floor, -4503599627370495.5L, -4503599627370496.0L); TEST_f_f (floor, -4503599627370496.25L, -4503599627370497.0L); TEST_f_f (floor, -4503599627370496.5L, -4503599627370497.0L); TEST_f_f (floor, -4503599627370496.75L, -4503599627370497.0L); TEST_f_f (floor, -4503599627370497.5L, -4503599627370498.0L); +# if LDBL_MANT_DIG > 100 + TEST_f_f (floor, -4503599627370494.5000000000001L, -4503599627370495.0L); + TEST_f_f (floor, -4503599627370495.5000000000001L, -4503599627370496.0L); + TEST_f_f (floor, -4503599627370496.5000000000001L, -4503599627370497.0L); +# endif TEST_f_f (floor, 9007199254740991.5L, 9007199254740991.0L); TEST_f_f (floor, 9007199254740992.25L, 9007199254740992.0L); @@ -2652,12 +2694,30 @@ floor_test (void) TEST_f_f (floor, 9007199254740992.75L, 9007199254740992.0L); TEST_f_f (floor, 9007199254740993.5L, 9007199254740993.0L); +# if LDBL_MANT_DIG > 100 + TEST_f_f (floor, 9007199254740991.0000000000001L, 9007199254740991.0L); + TEST_f_f (floor, 9007199254740992.0000000000001L, 9007199254740992.0L); + TEST_f_f (floor, 9007199254740993.0000000000001L, 9007199254740993.0L); + TEST_f_f (floor, 9007199254740991.5000000000001L, 9007199254740991.0L); + TEST_f_f (floor, 9007199254740992.5000000000001L, 9007199254740992.0L); + TEST_f_f (floor, 9007199254740993.5000000000001L, 9007199254740993.0L); +# endif + TEST_f_f (floor, -9007199254740991.5L, -9007199254740992.0L); TEST_f_f (floor, -9007199254740992.25L, -9007199254740993.0L); TEST_f_f (floor, -9007199254740992.5L, -9007199254740993.0L); TEST_f_f (floor, -9007199254740992.75L, -9007199254740993.0L); TEST_f_f (floor, -9007199254740993.5L, -9007199254740994.0L); +# if LDBL_MANT_DIG > 100 + TEST_f_f (floor, -9007199254740991.0000000000001L, -9007199254740992.0L); + TEST_f_f (floor, -9007199254740992.0000000000001L, -9007199254740993.0L); + TEST_f_f (floor, -9007199254740993.0000000000001L, -9007199254740994.0L); + TEST_f_f (floor, -9007199254740991.5000000000001L, -9007199254740992.0L); + TEST_f_f (floor, -9007199254740992.5000000000001L, -9007199254740993.0L); + TEST_f_f (floor, -9007199254740993.5000000000001L, -9007199254740994.0L); +# endif + TEST_f_f (floor, 72057594037927935.5L, 72057594037927935.0L); TEST_f_f (floor, 72057594037927936.25L, 72057594037927936.0L); TEST_f_f (floor, 72057594037927936.5L, 72057594037927936.0L); @@ -3971,6 +4031,12 @@ rint_test (void) TEST_f_f (rint, -2.5, -2.0); TEST_f_f (rint, -3.5, -4.0); TEST_f_f (rint, -4.5, -4.0); + TEST_f_f (rint, 0.1, 0.0); + TEST_f_f (rint, 0.25, 0.0); + TEST_f_f (rint, 0.625, 1.0); + TEST_f_f (rint, -0.1, -0.0); + TEST_f_f (rint, -0.25, -0.0); + TEST_f_f (rint, -0.625, -1.0); #ifdef TEST_LDOUBLE /* The result can only be represented in long double. */ TEST_f_f (rint, 4503599627370495.5L, 4503599627370496.0L); @@ -3978,12 +4044,38 @@ rint_test (void) TEST_f_f (rint, 4503599627370496.5L, 4503599627370496.0L); TEST_f_f (rint, 4503599627370496.75L, 4503599627370497.0L); TEST_f_f (rint, 4503599627370497.5L, 4503599627370498.0L); + +# if LDBL_MANT_DIG > 100 + TEST_f_f (rint, 4503599627370494.5000000000001L, 4503599627370495.0L); + TEST_f_f (rint, 4503599627370495.5000000000001L, 4503599627370496.0L); + TEST_f_f (rint, 4503599627370496.5000000000001L, 4503599627370497.0L); +# endif TEST_f_f (rint, -4503599627370495.5L, -4503599627370496.0L); TEST_f_f (rint, -4503599627370496.25L, -4503599627370496.0L); TEST_f_f (rint, -4503599627370496.5L, -4503599627370496.0L); TEST_f_f (rint, -4503599627370496.75L, -4503599627370497.0L); TEST_f_f (rint, -4503599627370497.5L, -4503599627370498.0L); + +# if LDBL_MANT_DIG > 100 + TEST_f_f (rint, -4503599627370494.5000000000001L, -4503599627370495.0L); + TEST_f_f (rint, -4503599627370495.5000000000001L, -4503599627370496.0L); + TEST_f_f (rint, -4503599627370496.5000000000001L, -4503599627370497.0L); + + TEST_f_f (rint, 9007199254740991.0000000000001L, 9007199254740991.0L); + TEST_f_f (rint, 9007199254740992.0000000000001L, 9007199254740992.0L); + TEST_f_f (rint, 9007199254740993.0000000000001L, 9007199254740993.0L); + TEST_f_f (rint, 9007199254740991.5000000000001L, 9007199254740992.0L); + TEST_f_f (rint, 9007199254740992.5000000000001L, 9007199254740993.0L); + TEST_f_f (rint, 9007199254740993.5000000000001L, 9007199254740994.0L); + + TEST_f_f (rint, -9007199254740991.0000000000001L, -9007199254740991.0L); + TEST_f_f (rint, -9007199254740992.0000000000001L, -9007199254740992.0L); + TEST_f_f (rint, -9007199254740993.0000000000001L, -9007199254740993.0L); + TEST_f_f (rint, -9007199254740991.5000000000001L, -9007199254740992.0L); + TEST_f_f (rint, -9007199254740992.5000000000001L, -9007199254740993.0L); + TEST_f_f (rint, -9007199254740993.5000000000001L, -9007199254740994.0L); +# endif TEST_f_f (rint, 9007199254740991.5L, 9007199254740992.0L); TEST_f_f (rint, 9007199254740992.25L, 9007199254740992.0L); @@ -4039,6 +4131,49 @@ rint_test_tonearest (void) TEST_f_f (rint, -1.0, -1.0); TEST_f_f (rint, -1.5, -2.0); TEST_f_f (rint, -2.0, -2.0); + TEST_f_f (rint, 0.1, 0.0); + TEST_f_f (rint, 0.25, 0.0); + TEST_f_f (rint, 0.625, 1.0); + TEST_f_f (rint, -0.1, -0.0); + TEST_f_f (rint, -0.25, -0.0); + TEST_f_f (rint, -0.625, -1.0); +#ifdef TEST_LDOUBLE + /* The result can only be represented in long double. */ + TEST_f_f (rint, 4503599627370495.5L, 4503599627370496.0L); + TEST_f_f (rint, 4503599627370496.25L, 4503599627370496.0L); + TEST_f_f (rint, 4503599627370496.5L, 4503599627370496.0L); + TEST_f_f (rint, 4503599627370496.75L, 4503599627370497.0L); + TEST_f_f (rint, 4503599627370497.5L, 4503599627370498.0L); +# if LDBL_MANT_DIG > 100 + TEST_f_f (rint, 4503599627370494.5000000000001L, 4503599627370495.0L); + TEST_f_f (rint, 4503599627370495.5000000000001L, 4503599627370496.0L); + TEST_f_f (rint, 4503599627370496.5000000000001L, 4503599627370497.0L); +# endif + TEST_f_f (rint, -4503599627370495.5L, -4503599627370496.0L); + TEST_f_f (rint, -4503599627370496.25L, -4503599627370496.0L); + TEST_f_f (rint, -4503599627370496.5L, -4503599627370496.0L); + TEST_f_f (rint, -4503599627370496.75L, -4503599627370497.0L); + TEST_f_f (rint, -4503599627370497.5L, -4503599627370498.0L); +# if LDBL_MANT_DIG > 100 + TEST_f_f (rint, -4503599627370494.5000000000001L, -4503599627370495.0L); + TEST_f_f (rint, -4503599627370495.5000000000001L, -4503599627370496.0L); + TEST_f_f (rint, -4503599627370496.5000000000001L, -4503599627370497.0L); + + TEST_f_f (rint, 9007199254740991.0000000000001L, 9007199254740991.0L); + TEST_f_f (rint, 9007199254740992.0000000000001L, 9007199254740992.0L); + TEST_f_f (rint, 9007199254740993.0000000000001L, 9007199254740993.0L); + TEST_f_f (rint, 9007199254740991.5000000000001L, 9007199254740992.0L); + TEST_f_f (rint, 9007199254740992.5000000000001L, 9007199254740993.0L); + TEST_f_f (rint, 9007199254740993.5000000000001L, 9007199254740994.0L); + + TEST_f_f (rint, -9007199254740991.0000000000001L, -9007199254740991.0L); + TEST_f_f (rint, -9007199254740992.0000000000001L, -9007199254740992.0L); + TEST_f_f (rint, -9007199254740993.0000000000001L, -9007199254740993.0L); + TEST_f_f (rint, -9007199254740991.5000000000001L, -9007199254740992.0L); + TEST_f_f (rint, -9007199254740992.5000000000001L, -9007199254740993.0L); + TEST_f_f (rint, -9007199254740993.5000000000001L, -9007199254740994.0L); +# endif +#endif } fesetround(save_round_mode); @@ -4066,6 +4201,49 @@ rint_test_towardzero (void) TEST_f_f (rint, -1.0, -1.0); TEST_f_f (rint, -1.5, -1.0); TEST_f_f (rint, -2.0, -2.0); + TEST_f_f (rint, 0.1, 0.0); + TEST_f_f (rint, 0.25, 0.0); + TEST_f_f (rint, 0.625, 0.0); + TEST_f_f (rint, -0.1, -0.0); + TEST_f_f (rint, -0.25, -0.0); + TEST_f_f (rint, -0.625, -0.0); +#ifdef TEST_LDOUBLE + /* The result can only be represented in long double. */ + TEST_f_f (rint, 4503599627370495.5L, 4503599627370495.0L); + TEST_f_f (rint, 4503599627370496.25L, 4503599627370496.0L); + TEST_f_f (rint, 4503599627370496.5L, 4503599627370496.0L); + TEST_f_f (rint, 4503599627370496.75L, 4503599627370496.0L); + TEST_f_f (rint, 4503599627370497.5L, 4503599627370497.0L); +# if LDBL_MANT_DIG > 100 + TEST_f_f (rint, 4503599627370494.5000000000001L, 4503599627370494.0L); + TEST_f_f (rint, 4503599627370495.5000000000001L, 4503599627370495.0L); + TEST_f_f (rint, 4503599627370496.5000000000001L, 4503599627370496.0L); +# endif + TEST_f_f (rint, -4503599627370495.5L, -4503599627370495.0L); + TEST_f_f (rint, -4503599627370496.25L, -4503599627370496.0L); + TEST_f_f (rint, -4503599627370496.5L, -4503599627370496.0L); + TEST_f_f (rint, -4503599627370496.75L, -4503599627370496.0L); + TEST_f_f (rint, -4503599627370497.5L, -4503599627370497.0L); +# if LDBL_MANT_DIG > 100 + TEST_f_f (rint, -4503599627370494.5000000000001L, -4503599627370494.0L); + TEST_f_f (rint, -4503599627370495.5000000000001L, -4503599627370495.0L); + TEST_f_f (rint, -4503599627370496.5000000000001L, -4503599627370496.0L); + + TEST_f_f (rint, 9007199254740991.0000000000001L, 9007199254740991.0L); + TEST_f_f (rint, 9007199254740992.0000000000001L, 9007199254740992.0L); + TEST_f_f (rint, 9007199254740993.0000000000001L, 9007199254740993.0L); + TEST_f_f (rint, 9007199254740991.5000000000001L, 9007199254740991.0L); + TEST_f_f (rint, 9007199254740992.5000000000001L, 9007199254740992.0L); + TEST_f_f (rint, 9007199254740993.5000000000001L, 9007199254740993.0L); + + TEST_f_f (rint, -9007199254740991.0000000000001L, -9007199254740991.0L); + TEST_f_f (rint, -9007199254740992.0000000000001L, -9007199254740992.0L); + TEST_f_f (rint, -9007199254740993.0000000000001L, -9007199254740993.0L); + TEST_f_f (rint, -9007199254740991.5000000000001L, -9007199254740991.0L); + TEST_f_f (rint, -9007199254740992.5000000000001L, -9007199254740992.0L); + TEST_f_f (rint, -9007199254740993.5000000000001L, -9007199254740993.0L); +# endif +#endif } fesetround(save_round_mode); @@ -4093,6 +4271,49 @@ rint_test_downward (void) TEST_f_f (rint, -1.0, -1.0); TEST_f_f (rint, -1.5, -2.0); TEST_f_f (rint, -2.0, -2.0); + TEST_f_f (rint, 0.1, 0.0); + TEST_f_f (rint, 0.25, 0.0); + TEST_f_f (rint, 0.625, 0.0); + TEST_f_f (rint, -0.1, -1.0); + TEST_f_f (rint, -0.25, -1.0); + TEST_f_f (rint, -0.625, -1.0); +#ifdef TEST_LDOUBLE + /* The result can only be represented in long double. */ + TEST_f_f (rint, 4503599627370495.5L, 4503599627370495.0L); + TEST_f_f (rint, 4503599627370496.25L, 4503599627370496.0L); + TEST_f_f (rint, 4503599627370496.5L, 4503599627370496.0L); + TEST_f_f (rint, 4503599627370496.75L, 4503599627370496.0L); + TEST_f_f (rint, 4503599627370497.5L, 4503599627370497.0L); +# if LDBL_MANT_DIG > 100 + TEST_f_f (rint, 4503599627370494.5000000000001L, 4503599627370494.0L); + TEST_f_f (rint, 4503599627370495.5000000000001L, 4503599627370495.0L); + TEST_f_f (rint, 4503599627370496.5000000000001L, 4503599627370496.0L); +# endif + TEST_f_f (rint, -4503599627370495.5L, -4503599627370496.0L); + TEST_f_f (rint, -4503599627370496.25L, -4503599627370497.0L); + TEST_f_f (rint, -4503599627370496.5L, -4503599627370497.0L); + TEST_f_f (rint, -4503599627370496.75L, -4503599627370497.0L); + TEST_f_f (rint, -4503599627370497.5L, -4503599627370498.0L); +# if LDBL_MANT_DIG > 100 + TEST_f_f (rint, -4503599627370494.5000000000001L, -4503599627370495.0L); + TEST_f_f (rint, -4503599627370495.5000000000001L, -4503599627370496.0L); + TEST_f_f (rint, -4503599627370496.5000000000001L, -4503599627370497.0L); + + TEST_f_f (rint, 9007199254740991.0000000000001L, 9007199254740991.0L); + TEST_f_f (rint, 9007199254740992.0000000000001L, 9007199254740992.0L); + TEST_f_f (rint, 9007199254740993.0000000000001L, 9007199254740993.0L); + TEST_f_f (rint, 9007199254740991.5000000000001L, 9007199254740991.0L); + TEST_f_f (rint, 9007199254740992.5000000000001L, 9007199254740992.0L); + TEST_f_f (rint, 9007199254740993.5000000000001L, 9007199254740993.0L); + + TEST_f_f (rint, -9007199254740991.0000000000001L, -9007199254740992.0L); + TEST_f_f (rint, -9007199254740992.0000000000001L, -9007199254740993.0L); + TEST_f_f (rint, -9007199254740993.0000000000001L, -9007199254740994.0L); + TEST_f_f (rint, -9007199254740991.5000000000001L, -9007199254740992.0L); + TEST_f_f (rint, -9007199254740992.5000000000001L, -9007199254740993.0L); + TEST_f_f (rint, -9007199254740993.5000000000001L, -9007199254740994.0L); +# endif +#endif } fesetround(save_round_mode); @@ -4120,6 +4341,49 @@ rint_test_upward (void) TEST_f_f (rint, -1.0, -1.0); TEST_f_f (rint, -1.5, -1.0); TEST_f_f (rint, -2.0, -2.0); + TEST_f_f (rint, 0.1, 1.0); + TEST_f_f (rint, 0.25, 1.0); + TEST_f_f (rint, 0.625, 1.0); + TEST_f_f (rint, -0.1, -0.0); + TEST_f_f (rint, -0.25, -0.0); + TEST_f_f (rint, -0.625, -0.0); +#ifdef TEST_LDOUBLE + /* The result can only be represented in long double. */ + TEST_f_f (rint, 4503599627370495.5L, 4503599627370496.0L); + TEST_f_f (rint, 4503599627370496.25L, 4503599627370497.0L); + TEST_f_f (rint, 4503599627370496.5L, 4503599627370497.0L); + TEST_f_f (rint, 4503599627370496.75L, 4503599627370497.0L); + TEST_f_f (rint, 4503599627370497.5L, 4503599627370498.0L); +# if LDBL_MANT_DIG > 100 + TEST_f_f (rint, 4503599627370494.5000000000001L, 4503599627370495.0L); + TEST_f_f (rint, 4503599627370495.5000000000001L, 4503599627370496.0L); + TEST_f_f (rint, 4503599627370496.5000000000001L, 4503599627370497.0L); +# endif + TEST_f_f (rint, -4503599627370495.5L, -4503599627370495.0L); + TEST_f_f (rint, -4503599627370496.25L, -4503599627370496.0L); + TEST_f_f (rint, -4503599627370496.5L, -4503599627370496.0L); + TEST_f_f (rint, -4503599627370496.75L, -4503599627370496.0L); + TEST_f_f (rint, -4503599627370497.5L, -4503599627370497.0L); +# if LDBL_MANT_DIG > 100 + TEST_f_f (rint, -4503599627370494.5000000000001L, -4503599627370494.0L); + TEST_f_f (rint, -4503599627370495.5000000000001L, -4503599627370495.0L); + TEST_f_f (rint, -4503599627370496.5000000000001L, -4503599627370496.0L); + + TEST_f_f (rint, 9007199254740991.0000000000001L, 9007199254740992.0L); + TEST_f_f (rint, 9007199254740992.0000000000001L, 9007199254740993.0L); + TEST_f_f (rint, 9007199254740993.0000000000001L, 9007199254740994.0L); + TEST_f_f (rint, 9007199254740991.5000000000001L, 9007199254740992.0L); + TEST_f_f (rint, 9007199254740992.5000000000001L, 9007199254740993.0L); + TEST_f_f (rint, 9007199254740993.5000000000001L, 9007199254740994.0L); + + TEST_f_f (rint, -9007199254740991.0000000000001L, -9007199254740991.0L); + TEST_f_f (rint, -9007199254740992.0000000000001L, -9007199254740992.0L); + TEST_f_f (rint, -9007199254740993.0000000000001L, -9007199254740993.0L); + TEST_f_f (rint, -9007199254740991.5000000000001L, -9007199254740991.0L); + TEST_f_f (rint, -9007199254740992.5000000000001L, -9007199254740992.0L); + TEST_f_f (rint, -9007199254740993.5000000000001L, -9007199254740993.0L); +# endif +#endif } fesetround(save_round_mode); @@ -4142,6 +4406,12 @@ round_test (void) TEST_f_f (round, -0.8L, -1.0); TEST_f_f (round, 1.5, 2.0); TEST_f_f (round, -1.5, -2.0); + TEST_f_f (round, 0.1, 0.0); + TEST_f_f (round, 0.25, 0.0); + TEST_f_f (round, 0.625, 1.0); + TEST_f_f (round, -0.1, -0.0); + TEST_f_f (round, -0.25, -0.0); + TEST_f_f (round, -0.625, -1.0); TEST_f_f (round, 2097152.5, 2097153); TEST_f_f (round, -2097152.5, -2097153); @@ -4151,13 +4421,23 @@ round_test (void) TEST_f_f (round, 4503599627370496.25L, 4503599627370496.0L); TEST_f_f (round, 4503599627370496.5L, 4503599627370497.0L); TEST_f_f (round, 4503599627370496.75L, 4503599627370497.0L); - TEST_f_f (round, 4503599627370497.5L, 4503599627370498.0L); + TEST_f_f (round, 4503599627370497.5L, 4503599627370498.0L); +# if LDBL_MANT_DIG > 100 + TEST_f_f (round, 4503599627370494.5000000000001L, 4503599627370495.0L); + TEST_f_f (round, 4503599627370495.5000000000001L, 4503599627370496.0L); + TEST_f_f (round, 4503599627370496.5000000000001L, 4503599627370497.0L); +# endif TEST_f_f (round, -4503599627370495.5L, -4503599627370496.0L); TEST_f_f (round, -4503599627370496.25L, -4503599627370496.0L); TEST_f_f (round, -4503599627370496.5L, -4503599627370497.0L); TEST_f_f (round, -4503599627370496.75L, -4503599627370497.0L); TEST_f_f (round, -4503599627370497.5L, -4503599627370498.0L); +# if LDBL_MANT_DIG > 100 + TEST_f_f (round, -4503599627370494.5000000000001L, -4503599627370495.0L); + TEST_f_f (round, -4503599627370495.5000000000001L, -4503599627370496.0L); + TEST_f_f (round, -4503599627370496.5000000000001L, -4503599627370497.0L); +# endif TEST_f_f (round, 9007199254740991.5L, 9007199254740992.0L); TEST_f_f (round, 9007199254740992.25L, 9007199254740992.0L); @@ -4171,6 +4451,22 @@ round_test (void) TEST_f_f (round, -9007199254740992.75L, -9007199254740993.0L); TEST_f_f (round, -9007199254740993.5L, -9007199254740994.0L); +# if LDBL_MANT_DIG > 100 + TEST_f_f (round, 9007199254740991.0000000000001L, 9007199254740991.0L); + TEST_f_f (round, 9007199254740992.0000000000001L, 9007199254740992.0L); + TEST_f_f (round, 9007199254740993.0000000000001L, 9007199254740993.0L); + TEST_f_f (round, 9007199254740991.5000000000001L, 9007199254740992.0L); + TEST_f_f (round, 9007199254740992.5000000000001L, 9007199254740993.0L); + TEST_f_f (round, 9007199254740993.5000000000001L, 9007199254740994.0L); + + TEST_f_f (round, -9007199254740991.0000000000001L, -9007199254740991.0L); + TEST_f_f (round, -9007199254740992.0000000000001L, -9007199254740992.0L); + TEST_f_f (round, -9007199254740993.0000000000001L, -9007199254740993.0L); + TEST_f_f (round, -9007199254740991.5000000000001L, -9007199254740992.0L); + TEST_f_f (round, -9007199254740992.5000000000001L, -9007199254740993.0L); + TEST_f_f (round, -9007199254740993.5000000000001L, -9007199254740994.0L); +# endif + TEST_f_f (round, 72057594037927935.5L, 72057594037927936.0L); TEST_f_f (round, 72057594037927936.25L, 72057594037927936.0L); TEST_f_f (round, 72057594037927936.5L, 72057594037927937.0L); @@ -4541,7 +4837,11 @@ trunc_test (void) TEST_f_f (trunc, 0, 0); TEST_f_f (trunc, minus_zero, minus_zero); + TEST_f_f (trunc, 0.1, 0); + TEST_f_f (trunc, 0.25, 0); TEST_f_f (trunc, 0.625, 0); + TEST_f_f (trunc, -0.1, minus_zero); + TEST_f_f (trunc, -0.25, minus_zero); TEST_f_f (trunc, -0.625, minus_zero); TEST_f_f (trunc, 1, 1); TEST_f_f (trunc, -1, -1); @@ -4565,11 +4865,23 @@ trunc_test (void) TEST_f_f (trunc, 4503599627370496.75L, 4503599627370496.0L); TEST_f_f (trunc, 4503599627370497.5L, 4503599627370497.0L); +# if LDBL_MANT_DIG > 100 + TEST_f_f (trunc, 4503599627370494.5000000000001L, 4503599627370494.0L); + TEST_f_f (trunc, 4503599627370495.5000000000001L, 4503599627370495.0L); + TEST_f_f (trunc, 4503599627370496.5000000000001L, 4503599627370496.0L); +# endif + TEST_f_f (trunc, -4503599627370495.5L, -4503599627370495.0L); TEST_f_f (trunc, -4503599627370496.25L, -4503599627370496.0L); TEST_f_f (trunc, -4503599627370496.5L, -4503599627370496.0L); TEST_f_f (trunc, -4503599627370496.75L, -4503599627370496.0L); TEST_f_f (trunc, -4503599627370497.5L, -4503599627370497.0L); + +# if LDBL_MANT_DIG > 100 + TEST_f_f (trunc, -4503599627370494.5000000000001L, -4503599627370494.0L); + TEST_f_f (trunc, -4503599627370495.5000000000001L, -4503599627370495.0L); + TEST_f_f (trunc, -4503599627370496.5000000000001L, -4503599627370496.0L); +# endif TEST_f_f (trunc, 9007199254740991.5L, 9007199254740991.0L); TEST_f_f (trunc, 9007199254740992.25L, 9007199254740992.0L); @@ -4577,12 +4889,30 @@ trunc_test (void) TEST_f_f (trunc, 9007199254740992.75L, 9007199254740992.0L); TEST_f_f (trunc, 9007199254740993.5L, 9007199254740993.0L); +# if LDBL_MANT_DIG > 100 + TEST_f_f (trunc, 9007199254740991.0000000000001L, 9007199254740991.0L); + TEST_f_f (trunc, 9007199254740992.0000000000001L, 9007199254740992.0L); + TEST_f_f (trunc, 9007199254740993.0000000000001L, 9007199254740993.0L); + TEST_f_f (trunc, 9007199254740991.5000000000001L, 9007199254740991.0L); + TEST_f_f (trunc, 9007199254740992.5000000000001L, 9007199254740992.0L); + TEST_f_f (trunc, 9007199254740993.5000000000001L, 9007199254740993.0L); +# endif + TEST_f_f (trunc, -9007199254740991.5L, -9007199254740991.0L); TEST_f_f (trunc, -9007199254740992.25L, -9007199254740992.0L); TEST_f_f (trunc, -9007199254740992.5L, -9007199254740992.0L); TEST_f_f (trunc, -9007199254740992.75L, -9007199254740992.0L); TEST_f_f (trunc, -9007199254740993.5L, -9007199254740993.0L); +# if LDBL_MANT_DIG > 100 + TEST_f_f (trunc, -9007199254740991.0000000000001L, -9007199254740991.0L); + TEST_f_f (trunc, -9007199254740992.0000000000001L, -9007199254740992.0L); + TEST_f_f (trunc, -9007199254740993.0000000000001L, -9007199254740993.0L); + TEST_f_f (trunc, -9007199254740991.5000000000001L, -9007199254740991.0L); + TEST_f_f (trunc, -9007199254740992.5000000000001L, -9007199254740992.0L); + TEST_f_f (trunc, -9007199254740993.5000000000001L, -9007199254740993.0L); +# endif + TEST_f_f (trunc, 72057594037927935.5L, 72057594037927935.0L); TEST_f_f (trunc, 72057594037927936.25L, 72057594037927936.0L); TEST_f_f (trunc, 72057594037927936.5L, 72057594037927936.0L); diff --git a/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c b/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c index 9666873..e99b0ba 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 416547c..8b1c976 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 e7e1b96..d055d65 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 a606548..035e4f5 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 2be5e28..4c4ae9b 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 51b1fea..1f4e33f 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 7fe84b8..0880e6e 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 14544af..d7bc47e 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 ae521fd..93663ad 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 7ae12a7..fd5fc0c 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 67518d0..a7efa3b 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 0000000..6cd6d0b --- /dev/null +++ b/sysdeps/powerpc/fpu/math_ldbl.h @@ -0,0 +1,189 @@ +#ifndef _MATH_PRIVATE_H_ +#error "Never use <math_ldbl.h> directly; include <math_private.h> instead." +#endif + +#include <sysdeps/ieee754/ldbl-128/math_ldbl.h> +#include <ieee754.h> + +static inline void +ldbl_extract_mantissa (int64_t *hi64, u_int64_t *lo64, int *exp, long double x) +{ + /* We have 105 bits of mantissa plus one implicit digit. Since + 106 bits are representable we use the first implicit digit for + the number before the decimal point and the second implicit bit + as bit 53 of the mantissa. */ + unsigned long long hi, lo; + int ediff; + union ibm_extended_long_double eldbl; + eldbl.d = x; + *exp = eldbl.ieee.exponent - IBM_EXTENDED_LONG_DOUBLE_BIAS; + + lo = ((long long)eldbl.ieee.mantissa2 << 32) | eldbl.ieee.mantissa3; + hi = ((long long)eldbl.ieee.mantissa0 << 32) | eldbl.ieee.mantissa1; + /* If the lower double is not a denomal or zero then set the hidden + 53rd bit. */ + if (eldbl.ieee.exponent2 > 0x001) + { + lo |= (1ULL << 52); + lo = lo << 7; /* pre-shift lo to match ieee854. */ + /* The lower double is normalized separately from the upper. We + may need to adjust the lower manitissa to reflect this. */ + ediff = eldbl.ieee.exponent - eldbl.ieee.exponent2; + if (ediff > 53) + lo = lo >> (ediff-53); + } + hi |= (1ULL << 52); + + if ((eldbl.ieee.negative != eldbl.ieee.negative2) + && ((eldbl.ieee.exponent2 != 0) && (lo != 0LL))) + { + hi--; + lo = (1ULL << 60) - lo; + if (hi < (1ULL << 52)) + { + /* we have a borrow from the hidden bit, so shift left 1. */ + hi = (hi << 1) | (lo >> 59); + lo = 0xfffffffffffffffLL & (lo << 1); + *exp = *exp - 1; + } + } + *lo64 = (hi << 60) | lo; + *hi64 = hi >> 4; +} + +static inline long double +ldbl_insert_mantissa (int sign, int exp, int64_t hi64, u_int64_t lo64) +{ + union ibm_extended_long_double u; + unsigned long hidden2, lzcount; + unsigned long long hi, lo; + + u.ieee.negative = sign; + u.ieee.negative2 = sign; + u.ieee.exponent = exp + IBM_EXTENDED_LONG_DOUBLE_BIAS; + u.ieee.exponent2 = exp-53 + IBM_EXTENDED_LONG_DOUBLE_BIAS; + /* Expect 113 bits (112 bits + hidden) right justified in two longs. + The low order 53 bits (52 + hidden) go into the lower double */ + lo = (lo64 >> 7)& ((1ULL << 53) - 1); + hidden2 = (lo64 >> 59) & 1ULL; + /* The high order 53 bits (52 + hidden) go into the upper double */ + hi = (lo64 >> 60) & ((1ULL << 11) - 1); + hi |= (hi64 << 4); + + if (lo != 0LL) + { + /* hidden2 bit of low double controls rounding of the high double. + If hidden2 is '1' then round up hi and adjust lo (2nd mantissa) + plus change the sign of the low double to compensate. */ + if (hidden2) + { + hi++; + u.ieee.negative2 = !sign; + lo = (1ULL << 53) - lo; + } + /* The hidden bit of the lo mantissa is zero so we need to + normalize the it for the low double. Shift it left until the + hidden bit is '1' then adjust the 2nd exponent accordingly. */ + + if (sizeof (lo) == sizeof (long)) + lzcount = __builtin_clzl (lo); + else if ((lo >> 32) != 0) + lzcount = __builtin_clzl ((long) (lo >> 32)); + else + lzcount = __builtin_clzl ((long) lo) + 32; + lzcount = lzcount - 11; + if (lzcount > 0) + { + int expnt2 = u.ieee.exponent2 - lzcount; + if (expnt2 >= 1) + { + /* Not denormal. Normalize and set low exponent. */ + lo = lo << lzcount; + u.ieee.exponent2 = expnt2; + } + else + { + /* Is denormal. */ + lo = lo << (lzcount + expnt2); + u.ieee.exponent2 = 0; + } + } + } + else + { + u.ieee.negative2 = 0; + u.ieee.exponent2 = 0; + } + + u.ieee.mantissa3 = lo & ((1ULL << 32) - 1); + u.ieee.mantissa2 = (lo >> 32) & ((1ULL << 20) - 1); + u.ieee.mantissa1 = hi & ((1ULL << 32) - 1); + u.ieee.mantissa0 = (hi >> 32) & ((1ULL << 20) - 1); + return u.d; +} + +/* gcc generates disgusting code to pack and unpack long doubles. + This tells gcc that pack/unpack is really a nop. We use fr1/fr2 + because those are the regs used to pass/return a single + long double arg. */ +static inline long double +ldbl_pack (double a, double aa) +{ + register long double x __asm__ ("fr1"); + register double xh __asm__ ("fr1"); + register double xl __asm__ ("fr2"); + xh = a; + xl = aa; + __asm__ ("" : "=f" (x) : "f" (xh), "f" (xl)); + return x; +} + +static inline void +ldbl_unpack (long double l, double *a, double *aa) +{ + register long double x __asm__ ("fr1"); + register double xh __asm__ ("fr1"); + register double xl __asm__ ("fr2"); + x = l; + __asm__ ("" : "=f" (xh), "=f" (xl) : "f" (x)); + *a = xh; + *aa = xl; +} + + +/* Convert a finite long double to canonical form. + Does not handle +/-Inf properly. */ +static inline void +ldbl_canonicalize (double *a, double *aa) +{ + double xh, xl; + + xh = *a + *aa; + xl = (*a - xh) + *aa; + *a = xh; + *aa = xl; +} + +/* Simple inline nearbyint (double) function . + Only works in the default rounding mode + but is useful in long double rounding functions. */ +static inline double +ldbl_nearbyint (double a) +{ + double two52 = 0x10000000000000LL; + + if (__builtin_expect ((__builtin_fabs (a) < two52), 1)) + { + if (__builtin_expect ((a > 0.0), 1)) + { + a += two52; + a -= two52; + } + else if (__builtin_expect ((a < 0.0), 1)) + { + a = two52 - a; + a = -(a - two52); + } + } + return a; +} diff --git a/sysdeps/powerpc/powerpc64/fpu/s_rintl.S b/sysdeps/powerpc/powerpc64/fpu/s_rintl.S deleted file mode 100644 index 2ca2d44..0000000 --- a/sysdeps/powerpc/powerpc64/fpu/s_rintl.S +++ /dev/null @@ -1,113 +0,0 @@ -/* Round to int long double floating-point values. - IBM extended format long double version. - Copyright (C) 2004, 2006 Free Software Foundation, Inc. - This file is part of the GNU C Library. - - The GNU C Library is free software; you can redistribute it and/or - modify it under the terms of the GNU Lesser General Public - License as published by the Free Software Foundation; either - version 2.1 of the License, or (at your option) any later version. - - The GNU C Library is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public - License along with the GNU C Library; if not, write to the Free - Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA - 02111-1307 USA. */ - -/* This has been coded in assembler because GCC makes such a mess of it - when it's coded in C. */ - -#include <sysdep.h> -#include <math_ldbl_opt.h> - - .section ".toc","aw" -.LC0: /* 2**52 */ - .tc FD_43300000_0[TC],0x4330000000000000 - .section ".text" - -ENTRY (__rintl) - lfd fp13,.LC0@toc(2) - fabs fp0,fp1 - fsub fp12,fp13,fp13 /* generate 0.0 */ - fabs fp9,fp2 - fcmpu cr7,fp0,fp13 /* if (fabs(x) > TWO52) */ - fcmpu cr6,fp1,fp12 /* if (x > 0.0) */ - bnl- cr7,.L2 - fmr fp2,fp12 - bng- cr6,.L1 - fadd fp1,fp1,fp13 /* x+= TWO52; */ - fsub fp1,fp1,fp13 /* x-= TWO52; */ - fabs fp1,fp1 /* if (x == 0.0) */ - blr /* x = 0.0; */ -.L1: - bnllr- cr6 /* if (x < 0.0) */ - fsub fp1,fp1,fp13 /* x-= TWO52; */ - fadd fp1,fp1,fp13 /* x+= TWO52; */ - fnabs fp1,fp1 /* if (x == 0.0) */ - blr /* x = -0.0; */ - -/* The high double is > TWO52 so we need to round the low double and - perhaps the high double. In this case we have to round the low - double and handle any adjustment to the high double that may be - caused by rounding (up). This is complicated by the fact that the - high double may already be rounded and the low double may have the - opposite sign to compensate.This gets a bit tricky so we use the - following algorithm: - - tau = floor(x_high/TWO52); - x0 = x_high - tau; - x1 = x_low + tau; - r1 = rint(x1); - y_high = x0 + r1; - y_low = x0 - y_high + r1; - return y; */ -.L2: - fcmpu cr7,fp9,fp13 /* if (|x_low| > TWO52) */ - fcmpu cr0,fp9,fp12 /* || (|x_low| == 0.0) */ - fcmpu cr5,fp2,fp12 /* if (x_low > 0.0) */ - bgelr- cr7 /* return x; */ - beqlr- cr0 - fdiv fp8,fp1,fp13 /* x_high/TWO52 */ - - bng- cr6,.L6 /* if (x > 0.0) */ - fctidz fp0,fp8 - fcfid fp8,fp0 /* tau = floor(x_high/TWO52); */ - fadd fp8,fp8,fp8 /* tau++; Make tau even */ - bng cr5,.L4 /* if (x_low > 0.0) */ - fmr fp3,fp1 - fmr fp4,fp2 - b .L5 -.L4: /* if (x_low < 0.0) */ - fsub fp3,fp1,fp8 /* x0 = x_high - tau; */ - fadd fp4,fp2,fp8 /* x1 = x_low + tau; */ -.L5: - fadd fp5,fp4,fp13 /* r1 = x1 + TWO52; */ - fsub fp5,fp5,fp13 /* r1 = r1 - TWO52; */ - b .L9 -.L6: /* if (x < 0.0) */ - fctidz fp0,fp8 - fcfid fp8,fp0 /* tau = floor(x_high/TWO52); */ - fadd fp8,fp8,fp8 /* tau++; Make tau even */ - bnl cr5,.L7 /* if (x_low < 0.0) */ - fmr fp3,fp1 - fmr fp4,fp2 - b .L8 -.L7: /* if (x_low > 0.0) */ - fsub fp3,fp1,fp8 /* x0 = x_high - tau; */ - fadd fp4,fp2,fp8 /* x1 = x_low + tau; */ -.L8: - fsub fp5,fp13,fp4 /* r1 = TWO52 - x1; */ - fsub fp0,fp5,fp13 /* r1 = - (r1 - TWO52); */ - fneg fp5,fp0 -.L9: - fadd fp1,fp3,fp5 /* y_high = x0 + r1; */ - fsub fp2,fp3,fp1 /* y_low = x0 - y_high + r1; */ - fadd fp2,fp2,fp5 - blr -END (__rintl) - -long_double_symbol (libm, __rintl, rintl) |