Fix ldbl-128ibm floorl for non-default rounding modes (bug 17899).

The ldbl-128ibm implementation of floorl is only correct in
round-to-nearest mode (in other modes, there are incorrect results and
overflow exceptions in some cases going beyond the incorrect signs of
zero results noted in bug 17899).  It is also unnecessarily
complicated, rounding both high and low parts to the nearest integer
and then adjusting for the semantics of floor, when it seems more
natural to take the floor of the high part (__floor optimized inline
versions can be used), and that of the low part if the high part is an
integer.  This patch makes it use that simpler approach, with a
canonicalization that works in all rounding modes (given that the only
way the result can be noncanonical is if taking the floor of a
negative noninteger low part increased its exponent).

Tested for powerpc, where over a thousand failures are removed from
test-ldouble.out (floorl problems affect many powl tests).

	[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.
This commit is contained in:
Joseph Myers 2016-02-18 21:31:10 +00:00
parent 656ee792a9
commit 1833769e19
3 changed files with 57 additions and 29 deletions

View File

@ -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.

View File

@ -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;
}
}
}

View File

@ -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);