Skip to content

Commit

Permalink
Fix ldbl-128ibm roundl for non-default rounding modes (bug 19594).
Browse files Browse the repository at this point in the history
The ldbl-128ibm implementation of roundl is only correct in
round-to-nearest mode (in other modes, there are incorrect results and
overflow exceptions in some cases).  This patch reimplements it along
the lines used for floorl, ceill and truncl, using __round on the high
part, and on the low part if the high part is an integer, and then
adjusting in the cases where this is incorrect.

Tested for powerpc.

	[BZ #19594]
	* sysdeps/ieee754/ldbl-128ibm/s_roundl.c (__roundl): Use __round
	on high and low parts then adjust result and use
	ldbl_canonicalize_int if needed.
  • Loading branch information
Joseph Myers committed Feb 18, 2016
1 parent e2310a2 commit b9a7633
Show file tree
Hide file tree
Showing 2 changed files with 39 additions and 36 deletions.
5 changes: 5 additions & 0 deletions ChangeLog
Original file line number Diff line number Diff line change
@@ -1,5 +1,10 @@
2016-02-18 Joseph Myers <joseph@codesourcery.com>

[BZ #19594]
* sysdeps/ieee754/ldbl-128ibm/s_roundl.c (__roundl): Use __round
on high and low parts then adjust result and use
ldbl_canonicalize_int if needed.

[BZ #19593]
* sysdeps/ieee754/ldbl-128ibm/s_truncl.c (__truncl): Use __trunc
on high part and __floor or __ceil on low part then use
Expand Down
70 changes: 34 additions & 36 deletions sysdeps/ieee754/ldbl-128ibm/s_roundl.c
Original file line number Diff line number Diff line change
Expand Up @@ -38,46 +38,44 @@ __roundl (long double x)
&& __builtin_isless (__builtin_fabs (xh),
__builtin_inf ()), 1))
{
double orig_xh;

/* Long double arithmetic, including the canonicalisation below,
only works in round-to-nearest mode. */

/* 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)
hi = __round (xh);
if (hi != xh)
{
if (xl > 0.0 || (xl == 0.0 && orig_xh > 0.0))
lo += 1.0;
/* The high part is not an integer; the low part only
affects the result if the high part is exactly half way
between two integers and the low part is nonzero with the
opposite sign. */
if (fabs (hi - xh) == 0.5)
{
if (xh > 0 && xl < 0)
xh = hi - 1;
else if (xh < 0 && xl > 0)
xh = hi + 1;
else
xh = hi;
}
else
xh = hi;
xl = 0;
}
else if (-xh == 0.5)
else
{
if (xl < 0.0 || (xl == 0.0 && orig_xh < 0.0))
lo -= 1.0;
/* The high part is a nonzero integer. */
lo = __round (xl);
if (fabs (lo - xl) == 0.5)
{
if (xh > 0 && xl < 0)
xl = lo + 1;
else if (xh < 0 && lo > 0)
xl = lo - 1;
else
xl = lo;
}
else
xl = lo;
xh = hi;
ldbl_canonicalize_int (&xh, &xl);
}

/* 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);
}

return ldbl_pack (xh, xl);
Expand Down

0 comments on commit b9a7633

Please sign in to comment.