User account creation filtered due to spam.

Bug 59666 - IBM long double arithmetic results invalid in non-default rounding modes
Summary: IBM long double arithmetic results invalid in non-default rounding modes
Status: UNCONFIRMED
Alias: None
Product: gcc
Classification: Unclassified
Component: target (show other bugs)
Version: 4.9.0
: P3 normal
Target Milestone: ---
Assignee: Not yet assigned to anyone
URL:
Keywords:
Depends on:
Blocks:
 
Reported: 2014-01-03 16:19 UTC by Joseph S. Myers
Modified: 2016-04-16 09:56 UTC (History)
3 users (show)

See Also:
Host:
Target: powerpc*-*-linux*
Build:
Known to work:
Known to fail:
Last reconfirmed:


Attachments

Note You need to log in before you can comment on or make changes to this bug.
Description Joseph S. Myers 2014-01-03 16:19:46 UTC
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.
Comment 1 Vincent Lefèvre 2016-03-01 13:53:12 UTC
(In reply to Joseph S. Myers from comment #0)
> 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 
[...]

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."

> 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 suppose that this is the case only when the FENV_ACCESS pragma is ON. ISO C99 §7.6.1 says (and ditto for C11): "If part of a program tests floating-point status flags, sets floating-point control modes, or runs under non-default mode settings, but was translated with the state for the FENV_ACCESS pragma ‘‘off’’, the behavior is undefined." The important point is "runs under non-default mode settings".

So, I suppose that any performance impact would be limited to programs that are translated with the FENV_ACCESS pragma set to ON. However this would still be an issue for libraries or users of these libraries, depending on what convention is chosen.
Comment 2 Vincent Lefèvre 2016-03-01 14:13:59 UTC
(In reply to Joseph S. Myers from comment #0)
> 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.

This would not be sufficient. Even without overflow/underflow, the rounding direction needs to be honored for the basic arithmetic operations.
Comment 3 joseph@codesourcery.com 2016-03-04 00:40:15 UTC
On Tue, 1 Mar 2016, vincent-gcc at vinc17 dot net wrote:

> > 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 suppose that this is the case only when the FENV_ACCESS pragma is ON. ISO C99

In GCC terms, that means -frounding-math.
Comment 4 joseph@codesourcery.com 2016-03-04 00:51:09 UTC
On Tue, 1 Mar 2016, vincent-gcc at vinc17 dot net wrote:

> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59666
> 
> --- Comment #2 from Vincent Lefèvre <vincent-gcc at vinc17 dot net> ---
> (In reply to Joseph S. Myers from comment #0)
> > 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.
> 
> This would not be sufficient. Even without overflow/underflow, the rounding
> direction needs to be honored for the basic arithmetic operations.

For a non-IEEE type where the accuracy is specified as within 3ulp 
(depending on the operation in question), there is no definition of what 
honouring the rounding mode means.  Such types fall within 5.2.4.2.2#6, 
"The accuracy of the floating-point operations (+, -, *, /) ... is 
implementation-defined ... The implementation may state that the accuracy 
is unknown.".

Note the proposed TC for DR#441: "The value of FLT_ROUNDS applies to all 
IEC 60559 types supported by the implementation, but need not apply to 
non-IEC 60559 types.".

Temporary setting of round-to-nearest and adjusting results suffices for 
uses in glibc; I have a libgcc patch to that effect (not suitable for 
inclusion because it doesn't condition things on a non-default configure 
option), together with an in-progress set of glibc patches to fix most of 
the glibc bugs shown up by test-ldouble failures.  I think roughly four 
such glibc fixes plus libm-test-ulps regeneration should allow 
test-ldouble to pass with patched libgcc (with unpatched libgcc needing 
also some extension of the existing set of XFAILs for this bug, before a 
default glibc build can pass test-ldouble on powerpc by default).  
(Passing with patched libgcc without XFAILs enables the use of random test 
generation to find bugs in glibc libm.  Passing with unpatched libgcc with 
XFAILs enabled the use of this test to detect regressions.)
Comment 5 Vincent Lefèvre 2016-03-04 02:35:51 UTC
(In reply to joseph@codesourcery.com from comment #4)
> For a non-IEEE type where the accuracy is specified as within 3ulp 
> (depending on the operation in question), there is no definition of what 
> honouring the rounding mode means.  Such types fall within 5.2.4.2.2#6, 
> "The accuracy of the floating-point operations (+, -, *, /) ... is 
> implementation-defined ... The implementation may state that the accuracy 
> is unknown.".

But directed rounding mode requirements are different from just accuracy requirements. For instance, with rounding toward negative infinity, one must have ROUND(x) ≤ x. This allows one to get guaranteed bounds and to do interval arithmetic, which is here more important than accuracy. If the rounding direction is not honored, one no longer gets guaranteed bounds/enclosures.

> Note the proposed TC for DR#441: "The value of FLT_ROUNDS applies to all 
> IEC 60559 types supported by the implementation, but need not apply to 
> non-IEC 60559 types.".

IMHO, this is a bad proposition. This would mean that a user who switches from double to long double could get incorrect results (enclosures) as a consequence.
Comment 6 joseph@codesourcery.com 2016-03-04 18:54:47 UTC
On Fri, 4 Mar 2016, vincent-gcc at vinc17 dot net wrote:

> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59666
> 
> --- Comment #5 from Vincent Lefèvre <vincent-gcc at vinc17 dot net> ---
> (In reply to joseph@codesourcery.com from comment #4)
> > For a non-IEEE type where the accuracy is specified as within 3ulp 
> > (depending on the operation in question), there is no definition of what 
> > honouring the rounding mode means.  Such types fall within 5.2.4.2.2#6, 
> > "The accuracy of the floating-point operations (+, -, *, /) ... is 
> > implementation-defined ... The implementation may state that the accuracy 
> > is unknown.".
> 
> But directed rounding mode requirements are different from just accuracy
> requirements. For instance, with rounding toward negative infinity, one must
> have ROUND(x) ≤ x. This allows one to get guaranteed bounds and to do interval

There are no such requirements in ISO C for non-IEEE types; you need to 
use IEEE types (such as float, double or __float128) for interval 
arithmetic.  Effectively, for non-IEEE types, arithmetic can be considered 
to work like most libm functions (where also there are no such 
requirements on how the results relate to the rounding mode).

The problems as far as I'm concerned with the IBM long double code in 
libgcc are where (even given options such as -frounding-math) it does 
worse than the default glibc conventions for most libm functions.
Comment 7 Vincent Lefèvre 2016-04-16 09:56:31 UTC
(In reply to joseph@codesourcery.com from comment #6)
> There are no such requirements in ISO C for non-IEEE types;

The standard says nowhere that the rounding directions that are provided by the implementation apply only to IEEE types. It even explicitly says that there may be additional implementation-defined rounding directions, thus not related to IEEE at all!