| Summary: | IBM long double arithmetic results invalid in non-default rounding modes | ||
|---|---|---|---|
| Product: | gcc | Reporter: | Joseph S. Myers <jsm28> |
| Component: | target | Assignee: | Not yet assigned to anyone <unassigned> |
| Status: | UNCONFIRMED --- | ||
| Severity: | normal | CC: | dje, egallager, iains, vincent-gcc, vital.had |
| Priority: | P3 | ||
| Version: | 4.9.0 | ||
| Target Milestone: | --- | ||
| Host: | Target: | powerpc*-*-linux* | |
| Build: | Known to work: | ||
| Known to fail: | Last reconfirmed: | ||
(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. (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. 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.
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.)
(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. 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.
(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! (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? (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... |
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.