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: 2014-06-02 22:30 UTC (History)
2 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.