diff options
Diffstat (limited to 'sysdeps/ieee754/ldbl-128ibm/s_llrintl.c')
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/s_llrintl.c | 168 |
1 files changed, 97 insertions, 71 deletions
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_llrintl.c b/sysdeps/ieee754/ldbl-128ibm/s_llrintl.c index bacf1bee0d..7228489098 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_llrintl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_llrintl.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,89 +37,114 @@ __llrintl (x) long double x; #endif { - static const double TWO52 = 4503599627370496.0L; - union ibm_extended_long_double u; - long long result = __LONG_LONG_MAX__; + double xh, xl; + long long res, hi, lo; + int save_round; - u.d = x; + ldbl_unpack (x, &xh, &xl); - if (fabs (u.dd[0]) < TWO52) + /* Limit the range of values handled by the conversion to long long. + We do this because we aren't sure whether that conversion properly + raises FE_INVALID. */ + if (__builtin_expect + ((__builtin_fabs (xh) <= -(double) (-__LONG_LONG_MAX__ - 1)), 1) +#if !defined (FE_INVALID) + || 1 +#endif + ) { - 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) + save_round = fegetround (); + + if (__builtin_expect ((xh == -(double) (-__LONG_LONG_MAX__ - 1)), 0)) { - high -= TWO52; - high += TWO52; - if (high == 0.0) high = -0.0; + /* When XH is 9223372036854775808.0, converting to long long will + overflow, resulting in an invalid operation. However, XL might + be negative and of sufficient magnitude that the overall long + double is in fact in range. Avoid raising an exception. In any + case we need to convert this value specially, because + the converted value is not exactly represented as a double + thus subtracting HI from XH suffers rounding error. */ + hi = __LONG_LONG_MAX__; + xh = 1.0; } - result = high; - } - else if (fabs (u.dd[1]) < TWO52 && u.dd[1] != 0.0) - { - double high, low, tau; - long long lowint; - /* 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) + 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. */ - - 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; + hi = (long long) xh; + xh -= hi; } - else if (u.dd[0] < 0.0) + ldbl_canonicalize (&xh, &xl); + + lo = (long long) xh; + + /* Peg at max/min values, assuming that the above conversions do so. + Strictly speaking, we can return anything for values that overflow, + but this is more useful. */ + res = hi + lo; + + /* This is just sign(hi) == sign(lo) && sign(res) != sign(hi). */ + if (__builtin_expect (((~(hi ^ lo) & (res ^ hi)) < 0), 0)) + goto overflow; + + xh -= lo; + ldbl_canonicalize (&xh, &xl); + + hi = res; + 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 (fabs (xh) < 0.5 + || (fabs (xh) == 0.5 + && ((xh > 0.0 && xl < 0.0) + || (xh < 0.0 && xl > 0.0) + || (xl == 0.0 && (res & 1) == 0)))) + return res; + + if (xh < 0.0) + res -= 1; + else + res += 1; + break; + + case FE_TOWARDZERO: + if (res > 0 && (xh < 0.0 || (xh == 0.0 && xl < 0.0))) + res -= 1; + else if (res < 0 && (xh > 0.0 || (xh == 0.0 && xl > 0.0))) + res += 1; + return res; + break; + + case FE_UPWARD: + if (xh > 0.0 || (xh == 0.0 && xl > 0.0)) + res += 1; + break; + + case FE_DOWNWARD: + if (xh < 0.0 || (xh == 0.0 && xl < 0.0)) + res -= 1; + break; } - lowint = low; - result = high; - result += lowint; + + if (__builtin_expect (((~(hi ^ (res - hi)) & (res ^ hi)) < 0), 0)) + goto overflow; + + return res; } else - result = u.dd[0]; + { + if (xh > 0.0) + hi = __LONG_LONG_MAX__; + else if (xh < 0.0) + hi = -__LONG_LONG_MAX__ - 1; + else + /* Nan */ + hi = 0; + } - return result; +overflow: +#ifdef FE_INVALID + feraiseexcept (FE_INVALID); +#endif + return hi; } long_double_symbol (libm, __llrintl, llrintl); |