diff options
| -rw-r--r-- | ChangeLog | 8 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/math_ldbl.h | 33 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/s_floorl.c | 45 |
3 files changed, 57 insertions, 29 deletions
diff --git a/ChangeLog b/ChangeLog index ffc8ea3661..6f79779923 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,11 @@ +2016-02-18 Joseph Myers <joseph@codesourcery.com> + + [BZ #17899] + * sysdeps/ieee754/ldbl-128ibm/math_ldbl.h (ldbl_canonicalize_int): + New function. + * sysdeps/ieee754/ldbl-128ibm/s_floorl.c (__floorl): Use __floor + on high and low parts then use ldbl_canonicalize_int if needed. + 2016-02-18 Adhemerval Zanella <adhemerval.zanella@linaro.org> * configure: Regenerated. diff --git a/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h b/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h index 051352f9f7..625ce00e13 100644 --- a/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h +++ b/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h @@ -230,3 +230,36 @@ ldbl_nearbyint (double a) } return a; } + +/* Canonicalize a result from an integer rounding function, in any + rounding mode. *A and *AA are finite and integers, with *A being + nonzero; if the result is not already canonical, *AA is plus or + minus a power of 2 that does not exceed the least set bit in + *A. */ +static inline void +ldbl_canonicalize_int (double *a, double *aa) +{ + int64_t ax, aax; + EXTRACT_WORDS64 (ax, *a); + EXTRACT_WORDS64 (aax, *aa); + int expdiff = ((ax >> 52) & 0x7ff) - ((aax >> 52) & 0x7ff); + if (expdiff <= 53) + { + if (expdiff == 53) + { + /* Half way between two double values; noncanonical iff the + low bit of A's mantissa is 1. */ + if ((ax & 1) != 0) + { + *a += 2 * *aa; + *aa = -*aa; + } + } + else + { + /* The sum can be represented in a single double. */ + *a += *aa; + *aa = 0; + } + } +} diff --git a/sysdeps/ieee754/ldbl-128ibm/s_floorl.c b/sysdeps/ieee754/ldbl-128ibm/s_floorl.c index 912230870a..a1469646b0 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_floorl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_floorl.c @@ -35,35 +35,22 @@ __floorl (long double x) && __builtin_isless (__builtin_fabs (xh), __builtin_inf ()), 1)) { - /* Long double arithmetic, including the canonicalisation below, - only works in round-to-nearest mode. */ - - /* 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); + hi = __floor (xh); + if (hi != xh) + { + /* The high part is not an integer; the low part does not + affect the result. */ + xh = hi; + xl = 0; + } + else + { + /* The high part is a nonzero integer. */ + lo = __floor (xl); + xh = hi; + xl = lo; + ldbl_canonicalize_int (&xh, &xl); + } } return ldbl_pack (xh, xl); |