[Bug c/61399] New: LDBL_MAX is incorrect with IBM long double format / overflow issues near large values

[vincent-gcc at vinc17 dot net]

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

            Bug ID: 61399
           Summary: LDBL_MAX is incorrect with IBM long double format /
                    overflow issues near large values
           Product: gcc
           Version: 4.7.2
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: c
          Assignee: unassigned at gcc dot gnu.org
          Reporter: vincent-gcc at vinc17 dot net

On PowerPC, which uses the IBM long double format for the type long double, one
gets:

LDBL_MANT_DIG = 106 (this is the precision of the FP model),
LDBL_MAX_EXP = 1024 (this is the maximum exponent of the FP model).

Even though the IBM long double format, a.k.a. double-double format, is not a
floating-point format, it is some superset of a floating-point format as
specified by the C standard (see C11, §5.2.4.2.2 from p1 to p3), where p = 106
and emax = 1024.

By definition, for radix b = 2, LDBL_MAX is the value (1 - 2^(-p)) * 2^emax
(see §5.2.4.2.2p12), which is the largest value representable in the
floating-point model.

However, with GCC 4.7.2 20121109 (Red Hat 4.7.2-8), I get:

LDBL_MAX = 0x1.fffffffffffff7ffffffffffff8p+1023

instead of 0x1.ffffffffffffffffffffffffff8p+1023.

The following program shows that this is not a display bug:

#include <stdio.h>
#include <float.h>

int main (void)
{
  long double dmax = DBL_MAX;
  printf ("%La\n", LDBL_MAX);
  printf ("%La\n", LDBL_MAX - dmax);
  printf ("%La\n", dmax + dmax * DBL_EPSILON / 4);
  printf ("%La\n", dmax + dmax * DBL_EPSILON / 2);
  return 0;
}

It outputs:

0x1.fffffffffffff7ffffffffffff8p+1023
0x1.ffffffffffffep+969
0x1.fffffffffffff7ffffffffffffcp+1023
inf

also showing that the arithmetic is buggy. I suppose that the high double is
the value rounded to double precision, but this rule is incorrect near the
largest values, due to the overflow.

One may choose to keep the behavior, i.e. consider that the high double is the
value rounded to double precision, but this means that the floating-point model
would need to be changed; otherwise some values are not representable, as shown
above.