From 5c68d401698a58cf7da150d9cce769fa6679ba5f Mon Sep 17 00:00:00 2001 From: Roland McGrath Date: Thu, 16 Mar 2006 11:47:24 +0000 Subject: [BZ #2423] 2006-03-07 Jakub Jelinek [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 Alan Modra * 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. --- sysdeps/ieee754/ldbl-128ibm/s_roundl.c | 124 ++++++++++++++------------------- 1 file changed, 51 insertions(+), 73 deletions(-) (limited to 'sysdeps/ieee754/ldbl-128ibm/s_roundl.c') diff --git a/sysdeps/ieee754/ldbl-128ibm/s_roundl.c b/sysdeps/ieee754/ldbl-128ibm/s_roundl.c index 7fe84b87de..0880e6ee22 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_roundl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_roundl.c @@ -22,7 +22,7 @@ when it's coded in C. */ #include -#include +#include #include #include #include @@ -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); -- cgit 1.4.1