* RFC: fixing ldbl-128ibm round-to-integer functions for non-default rounding modes
@ 2015-12-09 17:55 Joseph Myers
2016-01-21 17:27 ` Ping " Joseph Myers
0 siblings, 1 reply; 3+ messages in thread
From: Joseph Myers @ 2015-12-09 17:55 UTC (permalink / raw)
To: libc-alpha, munroesj
The ldbl-128ibm implementations of ceill, floorl, roundl and truncl
produce incorrect results and spurious overflow exceptions in
non-default rounding modes. In the case of floorl, over 1000
test-ldouble failures for other functions (mainly powl) arise from the
problems.
Unlike some such issues, these issues can readily be fixed without
involving any temporary changes of rounding mode (such modification of
floating-point state being known to involve performance problems).
The patch below for floorl illustrates the approach. The existing
function implementations are overly complicated, rounding high and low
parts to nearest integers and then adjusting for the semantics of the
function being implemented. But actually, for floor you can take the
floor of the high part - note that __floor has an optimized inline
version for newer processors that can be used. If the high part is
not an integer, the floor of the high part is the desired answer; if
the high part is an integer, you take the floor of the low part. In
the latter case you need to recanonicalize the result in case the low
part was negative and taking its floor increased its exponent, but
this can only produce a restricted set of noncanonical values that are
easy to canonicalize in a way that works in all rounding modes.
Any comments on this as an approach for fixing these issues for ceill,
floorl, roundl and truncl (all functions can be fixed in similar
ways)? Are there things that could be implemented in a different way
to improve performance? I think that with such fixes plus allowing
spurious underflow and inexact exceptions for ldbl-128ibm, and a real
fmal implementation, we should be able to get to the point where the
number of test failures arising from libgcc issues that can't be fixed
without significant performance cost is small enough that they can be
marked xfail-rounding:ldbl-128ibm to allow test-ldouble to pass.
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).
2015-12-09 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.
diff --git a/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h b/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h
index 051352f..625ce00 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 c911f19..c5973c6 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);
--
Joseph S. Myers
joseph@codesourcery.com
^ permalink raw reply [flat|nested] 3+ messages in thread
* Ping Re: RFC: fixing ldbl-128ibm round-to-integer functions for non-default rounding modes
2015-12-09 17:55 RFC: fixing ldbl-128ibm round-to-integer functions for non-default rounding modes Joseph Myers
@ 2016-01-21 17:27 ` Joseph Myers
2016-02-05 22:39 ` Joseph Myers
0 siblings, 1 reply; 3+ messages in thread
From: Joseph Myers @ 2016-01-21 17:27 UTC (permalink / raw)
To: libc-alpha, munroesj
Any comments on this
<https://sourceware.org/ml/libc-alpha/2015-12/msg00157.html> (for 2.24)?
--
Joseph S. Myers
joseph@codesourcery.com
^ permalink raw reply [flat|nested] 3+ messages in thread
* Re: Ping Re: RFC: fixing ldbl-128ibm round-to-integer functions for non-default rounding modes
2016-01-21 17:27 ` Ping " Joseph Myers
@ 2016-02-05 22:39 ` Joseph Myers
0 siblings, 0 replies; 3+ messages in thread
From: Joseph Myers @ 2016-02-05 22:39 UTC (permalink / raw)
To: libc-alpha, munroesj
On Thu, 21 Jan 2016, Joseph Myers wrote:
> Any comments on this
> <https://sourceware.org/ml/libc-alpha/2015-12/msg00157.html> (for 2.24)?
This is the patch fixing all of ceill, floorl, truncl and roundl that I
plan to apply for 2.24 (split by function and with appropriate ChangeLog
entries and bugs filed for ceill, truncl and roundl), given the lack of
comments.
diff --git a/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h b/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h
index 051352f..625ce00 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_ceill.c b/sysdeps/ieee754/ldbl-128ibm/s_ceill.c
index ac649b7..635fddc 100644
--- a/sysdeps/ieee754/ldbl-128ibm/s_ceill.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_ceill.c
@@ -35,42 +35,22 @@ __ceill (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 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);
+ hi = __ceil (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 = __ceil (xl);
+ xh = hi;
+ xl = lo;
+ ldbl_canonicalize_int (&xh, &xl);
+ }
}
return ldbl_pack (xh, xl);
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_floorl.c b/sysdeps/ieee754/ldbl-128ibm/s_floorl.c
index 9122308..a146964 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);
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_roundl.c b/sysdeps/ieee754/ldbl-128ibm/s_roundl.c
index 20813ed..b01510f 100644
--- a/sysdeps/ieee754/ldbl-128ibm/s_roundl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_roundl.c
@@ -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);
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_truncl.c b/sysdeps/ieee754/ldbl-128ibm/s_truncl.c
index df58b64..b7d4bb5 100644
--- a/sysdeps/ieee754/ldbl-128ibm/s_truncl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_truncl.c
@@ -35,50 +35,22 @@ __truncl (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 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 (orig_xh < 0.0)
+ hi = __trunc (xh);
+ if (hi != xh)
{
- if (xh > 0.0 || (xh == 0.0 && xl > 0.0))
- lo += 1.0;
+ /* The high part is not an integer; the low part does not
+ affect the result. */
+ xh = hi;
+ xl = 0;
}
else
{
- if (xh < 0.0 || (xh == 0.0 && xl < 0.0))
- lo -= 1.0;
+ /* The high part is a nonzero integer. */
+ lo = xh > 0 ? __floor (xl) : __ceil (xl);
+ xh = hi;
+ xl = lo;
+ 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);
-
- /* Ensure we return -0 rather than +0 when appropriate. */
- if (orig_xh < 0.0)
- xh = -__builtin_fabs (xh);
}
return ldbl_pack (xh, xl);
--
Joseph S. Myers
joseph@codesourcery.com
^ permalink raw reply [flat|nested] 3+ messages in thread
end of thread, other threads:[~2016-02-05 22:39 UTC | newest]
Thread overview: 3+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2015-12-09 17:55 RFC: fixing ldbl-128ibm round-to-integer functions for non-default rounding modes Joseph Myers
2016-01-21 17:27 ` Ping " Joseph Myers
2016-02-05 22:39 ` Joseph Myers
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).