diff options
Diffstat (limited to 'sysdeps/ieee754/ldbl-128ibm/s_ceill.c')
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/s_ceill.c | 125 |
1 files changed, 48 insertions, 77 deletions
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_ceill.c b/sysdeps/ieee754/ldbl-128ibm/s_ceill.c index a606548873..035e4f52ce 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_ceill.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_ceill.c @@ -19,7 +19,7 @@ 02111-1307 USA. */ #include <math.h> -#include <fenv.h> +#include <fenv_libc.h> #include <math_ldbl_opt.h> #include <float.h> #include <ieee754.h> @@ -34,87 +34,58 @@ __ceill (x) long double x; #endif { - static const double TWO52 = 4503599627370496.0L; - int mode = fegetround(); - union ibm_extended_long_double u; + double xh, xl, hi, lo; - u.d = x; + ldbl_unpack (x, &xh, &xl); - if (fabs (u.dd[0]) < TWO52) + /* Return Inf, Nan, +/-0 unchanged. */ + if (__builtin_expect (xh != 0.0 + && __builtin_isless (__builtin_fabs (xh), + __builtin_inf ()), 1)) { - double high = u.dd[0]; - fesetround(FE_UPWARD); - if (high > 0.0) - { - high += TWO52; - high -= TWO52; - if (high == -0.0) high = 0.0; - } - else if (high < 0.0) - { - high -= TWO52; - high += TWO52; - if (high == 0.0) high = -0.0; - } - u.dd[0] = high; - u.dd[1] = 0.0; - fesetround(mode); - } - else if (fabs (u.dd[1]) < TWO52 && u.dd[1] != 0.0) - { - double high, low; - /* In this case we have to round the low double and handle any - adjustment to the high double that may be caused by rounding - (up). This is complicated by the fact that the high double - may already be rounded and the low double may have the - opposite sign to compensate. */ - if (u.dd[0] > 0.0) - { - if (u.dd[1] > 0.0) - { - /* If the high/low doubles are the same sign then simply - round the low double. */ - high = u.dd[0]; - low = u.dd[1]; - } - else if (u.dd[1] < 0.0) - { - /* Else the high double is pre rounded and we need to - adjust for that. */ - high = nextafter (u.dd[0], 0.0); - low = u.dd[1] + (u.dd[0] - high); - } - fesetround(FE_UPWARD); - low += TWO52; - low -= TWO52; - fesetround(mode); - } - else if (u.dd[0] < 0.0) - { - if (u.dd[1] < 0.0) - { - /* If the high/low doubles are the same sign then simply - round the low double. */ - high = u.dd[0]; - low = u.dd[1]; - } - else if (u.dd[1] > 0.0) - { - /* Else the high double is pre rounded and we need to - adjust for that. */ - high = nextafter (u.dd[0], 0.0); - low = u.dd[1] + (u.dd[0] - high); - } - fesetround(FE_UPWARD); - low -= TWO52; - low += TWO52; - fesetround(mode); - } - u.dd[0] = high + low; - u.dd[1] = high - u.dd[0] + low; + double orig_xh; + int save_round = fegetround (); + + /* Long double arithmetic, including the canonicalisation below, + only works in round-to-nearest mode. */ + fesetround (FE_TONEAREST); + + /* Convert the high double to integer. */ + orig_xh = xh; + hi = ldbl_nearbyint (xh); + + /* Subtract integral high part from the value. */ + xh -= hi; + ldbl_canonicalize (&xh, &xl); + + /* Now convert the low double, adjusted for any remainder from the + high double. */ + lo = ldbl_nearbyint (xh); + + /* Adjust the result when the remainder is non-zero. nearbyint + rounds values to the nearest integer, and values halfway + between integers to the nearest even integer. ceill must + round towards +Inf. */ + xh -= lo; + ldbl_canonicalize (&xh, &xl); + + if (xh > 0.0 || (xh == 0.0 && xl > 0.0)) + lo += 1.0; + + /* Ensure the final value is canonical. In certain cases, + rounding causes hi,lo calculated so far to be non-canonical. */ + xh = hi; + xl = lo; + ldbl_canonicalize (&xh, &xl); + + /* Ensure we return -0 rather than +0 when appropriate. */ + if (orig_xh < 0.0) + xh = -__builtin_fabs (xh); + + fesetround (save_round); } - return u.d; + return ldbl_pack (xh, xl); } long_double_symbol (libm, __ceill, ceill); |