[Bug tree-optimization/107967] [13 regression] The gcc commit r13-3923 caused the glibc make check fails.

[jakub at gcc dot gnu.org]

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

--- Comment #8 from Jakub Jelinek <jakub at gcc dot gnu.org> ---
Actually, looking at:
288       // Be extra careful if there may be discrepancies between the
289       // compile and runtime results.
290       if ((mode_composite || (real_isneg (&inf) ? real_less (&result,
&value)
291                               : !real_less (&value, &result)))
292           && (inexact || !real_identical (&result, &value)))
I wonder if we didn't get it exactly wrong.
value is what is computed by real_arithmetic and result is that rounded to the
desired TYPE_MODE.
So, for non-mode_composite, if inf is negative, we want to find a low boundary.
Now, the value -> result conversion can round downward or upward or result can
be identical to value.
If real_less (&result, &value), it rounded it downward, if real_less (&value,
&result),
it rounded it upwards.
If we want low boundary and it was already rounded downward, we don't need to
round further.  So, I think we want
-         if ((mode_composite || (real_isneg (&inf) ? real_less (&result,
&value)
+         if ((mode_composite || (real_isneg (&inf) ? !real_less (&result,
&value)
                                  : !real_less (&value, &result)))
              && (inexact || !real_identical (&result, &value)))
which fixes the thinko/inconsistency.  Say on:
double
foo (void)
{
  const double huge = 1.0e+300;
  return huge * huge;
}

double
bar (void)
{
  const double huge = 1.0e+300;
  return huge * -huge;
}
current trunk optimizes foo to return Inf; but doesn't optimize bar to return
-Inf.

Now, at least for -fno-trapping-math -fno-rounding-math, I'd say we actually
want
return Inf; and return -Inf;, for -ftrapping-math per PR107608 resolution
compute
[Inf, Inf] or [-Inf, -Inf] ranges and only for -frounding-math use [DBL_MAX,
Inf]
and [-Inf, -DBL_MAX] ranges.
For that I guess we might want to special case
real_isinf (&result) && !real_isinf (&value) && !flag_rounding_math
cases or even just a subset of those when in round to even mode the computation
would
always round to infinity by far.

And another thing is mode_composite, if frange_nextafter is the right thing to
do
from infinity.  But I think (but can't find now where we say it) that double
double
is only supported in round to nearest mode, at least:
#include <fenv.h>
#include <stdio.h>

void
foo ()
{
  volatile double huge = 1.0e+308;
  volatile double inf = __builtin_inf ();
  fesetround (FE_DOWNWARD);
  volatile double r1 = huge + huge;
  volatile double r2 = huge * huge;
  volatile double r3 = huge + inf;
  volatile double r4 = r2 + huge;
  volatile double r5 = inf - 1.0;
  volatile double r6 = inf - huge;
  fesetround (FE_TONEAREST);
  printf ("%e %e %e %e %e %e\n", r1, r2, r3, r4, r5, r6);
}

void
bar ()
{
  volatile long double huge = 1.0e+308L;
  volatile long double inf = __builtin_infl ();
  fesetround (FE_DOWNWARD);
  volatile long double r1 = huge + huge;
  volatile long double r2 = huge * huge;
  volatile long double r3 = huge + inf;
  volatile long double r4 = r2 + huge;
  volatile long double r5 = inf - 1.0;
  volatile long double r6 = inf - huge;
  fesetround (FE_TONEAREST);
  printf ("%Le %Le %Le %Le %Le %Le\n", r1, r2, r3, r4, r5, r6);
}

int
main ()
{
  foo ();
  bar ();
  return 0;
}
goes wild on ppc64 and prints
1.797693e+308 1.797693e+308 inf 1.797693e+308 inf inf
1.797693e+308 -inf inf -inf inf inf
The -inf is quite far from +inf or LDBL_MAX.
Anyway, what I found is in libgcc/config/rs6000/ibm-ldouble-format is:
"Addition and subtraction are performed using library routines.  They
are not at present performed perfectly accurately, the result produced
will be within 1ulp of the range generated by adding or subtracting
1ulp from the input values, where a 'ulp' is 2^(e-106) given the
exponent 'e'.  In the presence of cancellation, this may be
arbitrarily inaccurate.  Subtraction is done by negation and addition.

Multiplication is also performed using a library routine.  Its result
will be within 2ulp of the correct result.

Division is also performed using a library routine.  Its result will
be within 3ulp of the correct result."
so for frange_arithmetics mode_composite we probably need to account for that
extra ulp on +-, 2ulps for * and 3ulps for /.