range-op-float: Fix up frange_arithmetic [PR107967]

  Hi!

The addition of PLUS/MINUS/MULT/RDIV_EXPR frange handlers causes
miscompilation of some of the libm routines, resulting in lots of
glibc test failures.  A part of them is purely PR107608 fold-overflow-1.c
etc. issues, say when the code does
  return -0.5 / 0.0;
and expects division by zero to be emitted, but we propagate -Inf
and avoid the operation.
But there are also various tests where we end up with different computed
value from the expected ones.  All those cases are like:
 is:          inf   inf
 should be:   1.18973149535723176502e+4932   0xf.fffffffffffffff0p+16380
 is:          inf   inf
 should be:   1.18973149535723176508575932662800701e+4932   0x1.ffffffffffffffffffffffffffffp+16383
 is:          inf   inf
 should be:   1.7976931348623157e+308   0x1.fffffffffffffp+1023
 is:          inf   inf
 should be:   3.40282346e+38   0x1.fffffep+127
and the corresponding source looks like:
static const double huge = 1.0e+300;
double whatever (...) {
...
  return huge * huge;
...
}
which for rounding to nearest or +inf should and does return +inf, but
for rounding to -inf or 0 should instead return nextafter (inf, -inf);
The rules IEEE754 has are that operations on +-Inf operands are exact
and produce +-Inf (except for the invalid ones that produce NaN) regardless
of rounding mode, while overflows:
"a) roundTiesToEven and roundTiesToAway carry all overflows to ∞ with the
sign of the intermediate result.
b) roundTowardZero carries all overflows to the format’s largest finite
number with the sign of the intermediate result.
c) roundTowardNegative carries positive overflows to the format’s largest
finite number, and carries negative overflows to −∞.
d) roundTowardPositive carries negative overflows to the format’s most
negative finite number, and carries positive overflows to +∞."

The behavior around overflows to -Inf or nextafter (-inf, inf) was actually
handled correctly, we'd construct [-INF, -MAX] ranges in those cases
because !real_less (&value, &result) in that case - value is finite
but larger in magnitude than what the format can represent (but GCC
internal's format can), while result is -INF in that case.
But for the overflows to +Inf or nextafter (inf, -inf) was handled
incorrectly, it tested real_less (&result, &value) rather than
!real_less (&result, &value), the former test is true when already the
rounding value -> result rounded down and in that case we shouldn't
round again, we should round down when it didn't.

So, in theory this could be fixed just by adding one ! character,
-  if ((mode_composite || (real_isneg (&inf) ? real_less (&result, &value)
+  if ((mode_composite || (real_isneg (&inf) ? !real_less (&result, &value)
 			  : !real_less (&value, &result)))
but the following patch goes further.  The distance between
nextafter (inf, -inf) and inf is large (infinite) and expressions like
1.0e+300 * 1.0e+300 always produce +inf in round to nearest mode by far,
so I think having low bound of nextafter (inf, -inf) in that case is
unnecessary.  But if it isn't multiplication but say addition and we are
inexact and very close to the boundary between rounding to nearest
maximum representable vs. rounding to nearest +inf, still using [MAX, +INF]
etc. ranges seems safer because we don't know exactly what we lost in the
inexact computation.

The following patch implements that.

Bootstrapped/regtested on x86_64-linux and i686-linux, ok for trunk?

2022-12-07  Jakub Jelinek  <jakub@redhat.com>

	PR tree-optimization/107967
	* range-op-float.cc (frange_arithmetic): Fix a thinko - if
	inf is negative, use nextafter if !real_less (&result, &value)
	rather than if real_less (&result, &value).  If result is +-INF
	while value is finite and -fno-rounding-math, don't do rounding
	if !inexact or if result is significantly above max representable
	value or below min representable value.

	* gcc.dg/pr107967-1.c: New test.
	* gcc.dg/pr107967-2.c: New test.
	* gcc.dg/pr107967-3.c: New test.


	Jakub
  

On 12/7/22 09:29, Jakub Jelinek wrote:
> Hi!
> 
> The addition of PLUS/MINUS/MULT/RDIV_EXPR frange handlers causes
> miscompilation of some of the libm routines, resulting in lots of
> glibc test failures.  A part of them is purely PR107608 fold-overflow-1.c
> etc. issues, say when the code does
>    return -0.5 / 0.0;
> and expects division by zero to be emitted, but we propagate -Inf
> and avoid the operation.
> But there are also various tests where we end up with different computed
> value from the expected ones.  All those cases are like:
>   is:          inf   inf
>   should be:   1.18973149535723176502e+4932   0xf.fffffffffffffff0p+16380
>   is:          inf   inf
>   should be:   1.18973149535723176508575932662800701e+4932   0x1.ffffffffffffffffffffffffffffp+16383
>   is:          inf   inf
>   should be:   1.7976931348623157e+308   0x1.fffffffffffffp+1023
>   is:          inf   inf
>   should be:   3.40282346e+38   0x1.fffffep+127
> and the corresponding source looks like:
> static const double huge = 1.0e+300;
> double whatever (...) {
> ...
>    return huge * huge;
> ...
> }
> which for rounding to nearest or +inf should and does return +inf, but
> for rounding to -inf or 0 should instead return nextafter (inf, -inf);
> The rules IEEE754 has are that operations on +-Inf operands are exact
> and produce +-Inf (except for the invalid ones that produce NaN) regardless
> of rounding mode, while overflows:
> "a) roundTiesToEven and roundTiesToAway carry all overflows to ∞ with the
> sign of the intermediate result.
> b) roundTowardZero carries all overflows to the format’s largest finite
> number with the sign of the intermediate result.
> c) roundTowardNegative carries positive overflows to the format’s largest
> finite number, and carries negative overflows to −∞.
> d) roundTowardPositive carries negative overflows to the format’s most
> negative finite number, and carries positive overflows to +∞."
> 
> The behavior around overflows to -Inf or nextafter (-inf, inf) was actually
> handled correctly, we'd construct [-INF, -MAX] ranges in those cases
> because !real_less (&value, &result) in that case - value is finite
> but larger in magnitude than what the format can represent (but GCC
> internal's format can), while result is -INF in that case.
> But for the overflows to +Inf or nextafter (inf, -inf) was handled
> incorrectly, it tested real_less (&result, &value) rather than
> !real_less (&result, &value), the former test is true when already the
> rounding value -> result rounded down and in that case we shouldn't
> round again, we should round down when it didn't.
> 
> So, in theory this could be fixed just by adding one ! character,
> -  if ((mode_composite || (real_isneg (&inf) ? real_less (&result, &value)
> +  if ((mode_composite || (real_isneg (&inf) ? !real_less (&result, &value)
>   			  : !real_less (&value, &result)))
> but the following patch goes further.  The distance between
> nextafter (inf, -inf) and inf is large (infinite) and expressions like
> 1.0e+300 * 1.0e+300 always produce +inf in round to nearest mode by far,
> so I think having low bound of nextafter (inf, -inf) in that case is
> unnecessary.  But if it isn't multiplication but say addition and we are
> inexact and very close to the boundary between rounding to nearest
> maximum representable vs. rounding to nearest +inf, still using [MAX, +INF]
> etc. ranges seems safer because we don't know exactly what we lost in the
> inexact computation.
> 
> The following patch implements that.
> 
> Bootstrapped/regtested on x86_64-linux and i686-linux, ok for trunk?
> 
> 2022-12-07  Jakub Jelinek  <jakub@redhat.com>
> 
> 	PR tree-optimization/107967
> 	* range-op-float.cc (frange_arithmetic): Fix a thinko - if
> 	inf is negative, use nextafter if !real_less (&result, &value)
> 	rather than if real_less (&result, &value).  If result is +-INF
> 	while value is finite and -fno-rounding-math, don't do rounding
> 	if !inexact or if result is significantly above max representable
> 	value or below min representable value.
> 
> 	* gcc.dg/pr107967-1.c: New test.
> 	* gcc.dg/pr107967-2.c: New test.
> 	* gcc.dg/pr107967-3.c: New test.
> 
> --- gcc/range-op-float.cc.jj	2022-12-06 10:25:16.594848892 +0100
> +++ gcc/range-op-float.cc	2022-12-06 20:53:47.751295689 +0100
> @@ -287,9 +287,64 @@ frange_arithmetic (enum tree_code code,
>   
>     // Be extra careful if there may be discrepancies between the
>     // compile and runtime results.
> -  if ((mode_composite || (real_isneg (&inf) ? real_less (&result, &value)
> -			  : !real_less (&value, &result)))
> -      && (inexact || !real_identical (&result, &value)))
> +  bool round = false;
> +  if (mode_composite)
> +    round = true;
> +  else if (real_isneg (&inf))
> +    {
> +      round = !real_less (&result, &value);
> +      if (real_isinf (&result, false)
> +	  && !real_isinf (&value)
> +	  && !flag_rounding_math)
> +	{
> +	  // Use just [+INF, +INF] rather than [MAX, +INF]
> +	  // even if value is larger than MAX and rounds to
> +	  // nearest to +INF.  Unless INEXACT is true, in
> +	  // that case we need some extra buffer.
> +	  if (!inexact)
> +	    round = false;
> +	  else
> +	    {
> +	      REAL_VALUE_TYPE tmp = result, tmp2;
> +	      frange_nextafter (mode, tmp, inf);
> +	      // TMP is at this point the maximum representable
> +	      // number.
> +	      real_arithmetic (&tmp2, MINUS_EXPR, &value, &tmp);
> +	      if (!real_isneg (&tmp2)
> +		  && (REAL_EXP (&tmp2) - REAL_EXP (&tmp)
> +		      >= 2 - REAL_MODE_FORMAT (mode)->p))
> +		round = false;
> +	    }
> +	}
> +    }

This chunk...

> +  else
> +    {
> +      round = !real_less (&value, &result);
> +      if (real_isinf (&result, true)
> +	  && !real_isinf (&value)
> +	  && !flag_rounding_math)
> +	{
> +	  // Use just [-INF, -INF] rather than [-INF, +MAX]
> +	  // even if value is smaller than -MAX and rounds to
> +	  // nearest to -INF.  Unless INEXACT is true, in
> +	  // that case we need some extra buffer.
> +	  if (!inexact)
> +	    round = false;
> +	  else
> +	    {
> +	      REAL_VALUE_TYPE tmp = result, tmp2;
> +	      frange_nextafter (mode, tmp, inf);
> +	      // TMP is at this point the minimum representable
> +	      // number.
> +	      real_arithmetic (&tmp2, MINUS_EXPR, &value, &tmp);
> +	      if (real_isneg (&tmp2)
> +		  && (REAL_EXP (&tmp2) - REAL_EXP (&tmp)
> +		      >= 2 - REAL_MODE_FORMAT (mode)->p))
> +		round = false;
> +	    }
> +	}
> +    }

...is quite similar to this one.  Could you abstract this?

Also:

 > +      if (real_isinf (&result, true)
 > +	  && !real_isinf (&value)
 > +	  && !flag_rounding_math)

Either abstract this out into a predicate since it's also repeated, or 
comment it.

Aldy

> +  if (round && (inexact || !real_identical (&result, &value)))
>       {
>         if (mode_composite)
>   	{
> --- gcc/testsuite/gcc.dg/pr107967-1.c.jj	2022-12-06 20:02:22.844086729 +0100
> +++ gcc/testsuite/gcc.dg/pr107967-1.c	2022-12-06 20:03:59.444683025 +0100
> @@ -0,0 +1,35 @@
> +/* PR tree-optimization/107967 */
> +/* { dg-do compile { target float64 } } */
> +/* { dg-options "-O2 -frounding-math -fno-trapping-math -fdump-tree-optimized" } */
> +/* { dg-add-options float64 } */
> +/* { dg-final { scan-tree-dump-not "return\[ \t]\*-?Inf;" "optimized" } } */
> +
> +_Float64
> +foo (void)
> +{
> +  const _Float64 huge = 1.0e+300f64;
> +  return huge * huge;
> +}
> +
> +_Float64
> +bar (void)
> +{
> +  const _Float64 huge = 1.0e+300f64;
> +  return huge * -huge;
> +}
> +
> +_Float64
> +baz (void)
> +{
> +  const _Float64 a = 0x1.fffffffffffffp+1023f64;
> +  const _Float64 b = 0x1.fffffffffffffp+970f64;
> +  return a + b;
> +}
> +
> +_Float64
> +qux (void)
> +{
> +  const _Float64 a = 0x1.fffffffffffffp+1023f64;
> +  const _Float64 b = 0x1.fffffffffffffp+969f64;
> +  return a + b;
> +}
> --- gcc/testsuite/gcc.dg/pr107967-2.c.jj	2022-12-06 20:02:29.683987331 +0100
> +++ gcc/testsuite/gcc.dg/pr107967-2.c	2022-12-06 20:03:48.685839355 +0100
> @@ -0,0 +1,35 @@
> +/* PR tree-optimization/107967 */
> +/* { dg-do compile { target float64 } } */
> +/* { dg-options "-O2 -fno-rounding-math -fno-trapping-math -fdump-tree-optimized" } */
> +/* { dg-add-options float64 } */
> +/* { dg-final { scan-tree-dump-times "return\[ \t]\*-?Inf;" 3 "optimized" } } */
> +
> +_Float64
> +foo (void)
> +{
> +  const _Float64 huge = 1.0e+300f64;
> +  return huge * huge;
> +}
> +
> +_Float64
> +bar (void)
> +{
> +  const _Float64 huge = 1.0e+300f64;
> +  return huge * -huge;
> +}
> +
> +_Float64
> +baz (void)
> +{
> +  const _Float64 a = 0x1.fffffffffffffp+1023f64;
> +  const _Float64 b = 0x1.fffffffffffffp+970f64;
> +  return a + b;
> +}
> +
> +_Float64
> +qux (void)
> +{
> +  const _Float64 a = 0x1.fffffffffffffp+1023f64;
> +  const _Float64 b = 0x1.fffffffffffffp+969f64;
> +  return a + b;
> +}
> --- gcc/testsuite/gcc.dg/pr107967-3.c.jj	2022-12-06 20:29:35.243370388 +0100
> +++ gcc/testsuite/gcc.dg/pr107967-3.c	2022-12-06 20:53:16.553748313 +0100
> @@ -0,0 +1,53 @@
> +/* PR tree-optimization/107967 */
> +/* { dg-do compile { target float64 } } */
> +/* { dg-options "-O2 -fno-rounding-math -fno-trapping-math -fdump-tree-optimized" } */
> +/* { dg-add-options float64 } */
> +/* { dg-final { scan-tree-dump-times "return\[ \t]\*-?Inf;" 3 "optimized" } } */
> +
> +_Float64
> +foo (_Float64 x)
> +{
> +  if (x >= 1.0e+300f64)
> +    ;
> +  else
> +    __builtin_unreachable ();
> +  return x * x;
> +}
> +
> +_Float64
> +bar (_Float64 x)
> +{
> +  if (x >= 1.0e+300f64)
> +    ;
> +  else
> +    __builtin_unreachable ();
> +  return x * -x;
> +}
> +
> +_Float64
> +baz (_Float64 a, _Float64 b)
> +{
> +  if (a >= 0x1.fffffffffffffp+1023f64)
> +    ;
> +  else
> +    __builtin_unreachable ();
> +  if (b >= 0x1.p+972f64)
> +    ;
> +  else
> +    __builtin_unreachable ();
> +  return a + b;
> +}
> +
> +_Float64
> +qux (_Float64 a, _Float64 b)
> +{
> +  if (a >= 0x1.fffffffffffffp+1023f64)
> +    ;
> +  else
> +    __builtin_unreachable ();
> +  if (b >= 0x1.fffffffffffffp+969f64)
> +    ;
> +  else
> +    __builtin_unreachable ();
> +  return a + b;
> +}
> 
> 	Jakub
>