[Bug target/59666] New: IBM long double arithmetic results invalid in non-default rounding modes

["jsm28 at gcc dot gnu.org" <gcc-bugzilla@gcc.gnu.org>]

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.