http://gcc.gnu.org/bugzilla/show_bug.cgi?id=59666
Bug ID: 59666
Summary: IBM long double arithmetic results invalid in
non-default rounding modes
Product: gcc
Version: 4.9.0
Status: UNCONFIRMED
Severity: normal
Priority: P3
Component: target
Assignee: unassigned at gcc dot gnu.org
Reporter: jsm28 at gcc dot gnu.org
CC: dje at gcc dot gnu.org
Target: powerpc*-*-linux*
The IBM long double functions in ibm-ldouble.c, when called in rounding modes
other than round-to-nearest, can produce results that are invalid (do not
satisfy the requirements in ibm-ldouble-format on what pairs of doubles make a
valid long double value, in particular as regards the high part being equal to
the sum of the two parts rounded to nearest). For example:
#include <fenv.h>
#include <float.h>
#include <stdio.h>
volatile long double a = LDBL_MAX;
int
main (void)
{
fesetround (FE_TOWARDZERO);
union u { long double ld; double d[2]; } x;
volatile long double r = a * a;
x.ld = a;
printf ("LDBL_MAX: %a %a\n", x.d[0], x.d[1]);
x.ld = r;
printf ("LDBL_MAX * LDBL_MAX: %a %a\n", x.d[0], x.d[1]);
return 0;
}
prints
LDBL_MAX: 0x1.fffffffffffffp+1023 0x1.ffffffffffffep+969
LDBL_MAX * LDBL_MAX: 0x1.fffffffffffffp+1023 0x1.fffffffffffffp+1023
where the value for LDBL_MAX * LDBL_MAX is not a valid long double at all.
(This isn't limited to overflow cases, although they may produce the greatest
errors; e.g. 0x1.fffffffffffffp-1L * 0x1.fffffffffffffp-1L in FE_UPWARD mode
produces (0x1.fffffffffffffp-1, -0x1.fffffffffffffp-54), where the high part is
not the sum of the two parts rounded to nearest.)
ISO C does not allow for arithmetic operations simply not working - producing
invalid results - for some types and rounding modes, although for non-IEEE
types they need not be correctly rounding.
I think the right approach for a fix will probably involve setting
round-to-nearest temporarily within the functions, then adjusting overflowing
and underflowing results based on the original rounding mode. I don't know
what performance impact that might have, and whether it might be best to avoid
that performance impact in the common case by having separate __gcc_*_round
functions that deal with saving / restoring the rounding mode and are only used
if -frounding-math, with the existing functions (not handling rounding modes)
being used otherwise.
