[Bug target/25661] New: Wrong long double to float conversion

[Bug target/25661] New: Wrong long double to float conversion

[amodra at bigpond dot net dot au]

Conversion from 128-bit long double to float goes via double.  This means the
conversions suffers from double rounding, and incorrect results may be
calculated for long double values near the midpoint of two adjacent float
values.

For example, f and g should have different values in the following.

#ifndef DO_DOUBLE
#define FP float
#define FP_MANT 24
#else
#define FP double
#define FP_MANT 53
#endif

volatile long double x;
volatile long double y;

int main (void)
{
  FP f, g;
  x = (1ull << 63) + (1ull << (63 - FP_MANT));
  y = 1.0;
  f = x + y;
  g = x - y;
  return f == g;
}


-- 
           Summary: Wrong long double to float conversion
           Product: gcc
           Version: 4.2.0
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: target
        AssignedTo: unassigned at gcc dot gnu dot org
        ReportedBy: amodra at bigpond dot net dot au
GCC target triplet: powerpc64-*-*


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=25661

[Bug target/25661] Wrong long double to float conversion

[pinskia at gcc dot gnu dot org]

------- Comment #1 from pinskia at gcc dot gnu dot org  2006-01-04 15:17 -------
Hmm, actually there is no rounding from TF to DF except for the fact that upper
potion of the TF mode is taken only.

I am thinking this testcase is invalid for 128bit IBM long double format.

IBM 128bit long double format is not IEEE 128bit format at all.

it means that a+b.


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=25661

[Bug target/25661] Wrong long double to float conversion

[joseph at codesourcery dot com]

------- Comment #2 from joseph at codesourcery dot com  2006-01-04 16:47 -------
Subject: Re:  Wrong long double to float conversion

On Wed, 4 Jan 2006, pinskia at gcc dot gnu dot org wrote:

> Hmm, actually there is no rounding from TF to DF except for the fact that upper
> potion of the TF mode is taken only.
> 
> I am thinking this testcase is invalid for 128bit IBM long double format.
> 
> IBM 128bit long double format is not IEEE 128bit format at all.

C99 Annex F allows non-IEC 60559 formats for long double and it seems 
conversions from such types to float must follow the IEC 60559 rules.  Do 
you have a pointer to a specification for double+double long double which 
does something other than compute to infinite precision then round to 
nearest representable value according to current rounding mode?


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=25661

[Bug target/25661] Wrong long double to float conversion

[pinskia at gcc dot gnu dot org]

------- Comment #3 from pinskia at gcc dot gnu dot org  2006-01-04 17:04 -------
http://publib16.boulder.ibm.com/pseries/en_US/aixprggd/genprogc/128bit_long_double_floating-point_datatype.htm


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=25661

[Bug target/25661] Wrong long double to float conversion

[pinskia at gcc dot gnu dot org]

------- Comment #4 from pinskia at gcc dot gnu dot org  2006-01-04 17:13 -------
>From that page:
The 128-bit implementation differs from the IEEE standard for long double in
the following ways:

Supports only round-to-nearest mode. If the application changes the rounding
mode, results are undefined.
Does not fully support the IEEE special numbers NaN and INF.
Does not support IEEE status flags for overflow, underflow, and other
conditions. These flags have no meaning for the 128-bit long double
inplementation.


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=25661

[Bug target/25661] Wrong long double to float conversion

[dje at gcc dot gnu dot org]

------- Comment #5 from dje at gcc dot gnu dot org  2006-01-05 18:26 -------
Andrew, you are confusing the issue by quoting the AIX documentation.  Darwin
and Linux implement different semantics from AIX.  On Darwin and Linux, the
most significant doubleword always is correctly rounded, so it contains the
properly rounded double value that is rounded to single precision.  On AIX, the
properly rounded double precision value is computed by summing the two
doublewords.


-- 

dje at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |dje at gcc dot gnu dot org


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=25661

[Bug target/25661] Wrong long double to float conversion

[amodra at bigpond dot net dot au]

------- Comment #6 from amodra at bigpond dot net dot au  2006-01-06 09:31 -------
In answer to comment #1, using just the upper double of the TFmode long double
*is* rounding to double.  As David says in #5, the upper double is the value of
the long double rounded to double..


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=25661

[Bug target/25661] Wrong long double to float conversion

[pinskia at gcc dot gnu dot org]

------- Comment #7 from pinskia at gcc dot gnu dot org  2006-01-27 18:03 -------
But this is by design and not a bug in GCC.


-- 

pinskia at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|UNCONFIRMED                 |RESOLVED
         Resolution|                            |INVALID


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=25661

[Bug target/25661] Wrong long double to float conversion

[amodra at bigpond dot net dot au]

------- Comment #8 from amodra at bigpond dot net dot au  2006-01-28 00:17 -------
There is no reason to say that wrong long double -> float conversion code
emitted by gcc is some fault of the long double design.  It is relatively easy
to convert correctly:

float long_double_to_float (long double a)
{
  union { long double ld; double d[2]; } u;
  double xh, xl, tmp;
  float res;

  u.ld = a;
  xh = u.d[0];
  xl = u.d[1];

  /* Convert high double.  */
  res = (float) xh;
  /* Subtract from long double, and canonicalise result.  */
  xh -= res;
  tmp = xh + xl;
  xl = (xh - tmp) + xl;
  xh = tmp;
  /* Convert remainder.  */
  res += (float) xh;
  return res;
}

There are some obvious simplifications to be made to the above function.  I
left the dead code there to show a technique that is useful in many long double
algorithms, ie. operate on high double, subtract result and canonicalise,
operate on low double (which has moved to xh by virtue of canonicalisation).


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=25661