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

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

[jsm28 at gcc dot 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.

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

[vital.had at gmail dot com]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59666

Sergey Fedorov <vital.had at gmail dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |vital.had at gmail dot com

--- Comment #8 from Sergey Fedorov <vital.had at gmail dot com> ---
(In reply to Vincent Lefèvre from comment #1)
> (In reply to Joseph S. Myers from comment #0)
> It seems to be like that "by design" (though this is not satisfactory) and
> part of the ppc64 ABI for instance:
> 
> http://refspecs.linuxfoundation.org/ELF/ppc64/PPC-elf64abi.html#PREC
> 
> "The software support is restricted to round-to-nearest mode. Programs that
> use extended precision must ensure that this rounding mode is in effect when
> extended-precision calculations are performed."

Also true for AIX:
https://www.ibm.com/docs/sr/xcafbg/9.0.0?topic=SS3KZ4_9.0.0/com.ibm.xlf111.bg.doc/xlfopg/fp-overview.html

Does anyone know whether this is also true for Darwin on PowerPC though?

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

[egallager at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59666

Eric Gallager <egallager at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |egallager at gcc dot gnu.org,
                   |                            |iains at gcc dot gnu.org

--- Comment #9 from Eric Gallager <egallager at gcc dot gnu.org> ---
(In reply to Sergey Fedorov from comment #8)
> (In reply to Vincent Lefèvre from comment #1)
> > (In reply to Joseph S. Myers from comment #0)
> > It seems to be like that "by design" (though this is not satisfactory) and
> > part of the ppc64 ABI for instance:
> > 
> > http://refspecs.linuxfoundation.org/ELF/ppc64/PPC-elf64abi.html#PREC
> > 
> > "The software support is restricted to round-to-nearest mode. Programs that
> > use extended precision must ensure that this rounding mode is in effect when
> > extended-precision calculations are performed."
> 
> Also true for AIX:
> https://www.ibm.com/docs/sr/xcafbg/9.0.0?topic=SS3KZ4_9.0.0/com.ibm.xlf111.
> bg.doc/xlfopg/fp-overview.html
> 
> Does anyone know whether this is also true for Darwin on PowerPC though?

I don't, but Iain might...