range-op: Implement floating point multiplication fold_range [PR107569]

Hi!

The following patch implements frange multiplication, including the
special case of x * x.  The callers don't tell us that it is x * x,
just that it is either z = x * x or if (x == y) z = x * y;
For irange that makes no difference, but for frange it can mean
x is -0.0 and y is 0.0 if they have the same range that includes both
signed and unsigned zeros, so we need to assume result could be -0.0.

The patch causes one regression:
+FAIL: gcc.dg/fold-overflow-1.c scan-assembler-times 2139095040 2
but that is already tracked in PR107608 and affects not just the newly
added multiplication, but addition and other floating point operations
(and doesn't seem like a ranger bug but dce or whatever else).

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

Once we have division and the reverse ops for all of these, perhaps
we can do some cleanups to share common code, but the way I have division
now partly written doesn't show up many commonalities.  Multiplication
is simple, division is a nightmare.

2022-11-10  Jakub Jelinek  <jakub@redhat.com>

	PR tree-optimization/107569
	PR tree-optimization/107591
	* range-op.h (range_operator_float::rv_fold): Add relation_kind
	argument.
	* range-op-float.cc (range_operator_float::fold_range): Name
	last argument trio and pass trio.op1_op2 () as last argument to
	rv_fold.
	(range_operator_float::rv_fold): Add relation_kind argument.
	(foperator_plus::rv_fold, foperator_minus::rv_fold): Likewise.
	(frange_mult): New function.
	(foperator_mult): New class.
	(floating_op_table::floating_op_table): Use foperator_mult for
	MULT_EXPR.


	Jakub

On 11/10/22 14:44, Jakub Jelinek wrote:
> Hi!
> 
> The following patch implements frange multiplication, including the
> special case of x * x.  The callers don't tell us that it is x * x,
> just that it is either z = x * x or if (x == y) z = x * y;
> For irange that makes no difference, but for frange it can mean
> x is -0.0 and y is 0.0 if they have the same range that includes both
> signed and unsigned zeros, so we need to assume result could be -0.0.
> 
> The patch causes one regression:
> +FAIL: gcc.dg/fold-overflow-1.c scan-assembler-times 2139095040 2
> but that is already tracked in PR107608 and affects not just the newly
> added multiplication, but addition and other floating point operations
> (and doesn't seem like a ranger bug but dce or whatever else).
> 
> Bootstrapped/regtested on x86_64-linux and i686-linux, ok for trunk?
> 
> Once we have division and the reverse ops for all of these, perhaps
> we can do some cleanups to share common code, but the way I have division
> now partly written doesn't show up many commonalities.  Multiplication
> is simple, division is a nightmare.

Thanks for tackling this.  I'm happy you think multiplication is simple, 
cause all the cross product operators make my head spin.

> 
> 2022-11-10  Jakub Jelinek  <jakub@redhat.com>
> 
> 	PR tree-optimization/107569
> 	PR tree-optimization/107591
> 	* range-op.h (range_operator_float::rv_fold): Add relation_kind
> 	argument.
> 	* range-op-float.cc (range_operator_float::fold_range): Name
> 	last argument trio and pass trio.op1_op2 () as last argument to
> 	rv_fold.
> 	(range_operator_float::rv_fold): Add relation_kind argument.
> 	(foperator_plus::rv_fold, foperator_minus::rv_fold): Likewise.
> 	(frange_mult): New function.
> 	(foperator_mult): New class.
> 	(floating_op_table::floating_op_table): Use foperator_mult for
> 	MULT_EXPR.
> 
> --- gcc/range-op.h.jj	2022-11-10 00:55:09.430219763 +0100
> +++ gcc/range-op.h	2022-11-10 11:30:33.594114939 +0100
> @@ -123,7 +123,8 @@ public:
>   			const REAL_VALUE_TYPE &lh_lb,
>   			const REAL_VALUE_TYPE &lh_ub,
>   			const REAL_VALUE_TYPE &rh_lb,
> -			const REAL_VALUE_TYPE &rh_ub) const;
> +			const REAL_VALUE_TYPE &rh_ub,
> +			relation_kind) const;
>     // Unary operations have the range of the LHS as op2.
>     virtual bool fold_range (irange &r, tree type,
>   			   const frange &lh,
> --- gcc/range-op-float.cc.jj	2022-11-10 00:55:09.318221259 +0100
> +++ gcc/range-op-float.cc	2022-11-10 11:31:29.040359082 +0100
> @@ -51,7 +51,7 @@ along with GCC; see the file COPYING3.
>   bool
>   range_operator_float::fold_range (frange &r, tree type,
>   				  const frange &op1, const frange &op2,
> -				  relation_trio) const
> +				  relation_trio trio) const
>   {
>     if (empty_range_varying (r, type, op1, op2))
>       return true;
> @@ -65,7 +65,7 @@ range_operator_float::fold_range (frange
>     bool maybe_nan;
>     rv_fold (lb, ub, maybe_nan, type,
>   	   op1.lower_bound (), op1.upper_bound (),
> -	   op2.lower_bound (), op2.upper_bound ());
> +	   op2.lower_bound (), op2.upper_bound (), trio.op1_op2 ());
>   
>     // Handle possible NANs by saturating to the appropriate INF if only
>     // one end is a NAN.  If both ends are a NAN, just return a NAN.
> @@ -103,8 +103,8 @@ range_operator_float::rv_fold (REAL_VALU
>   			       const REAL_VALUE_TYPE &lh_lb ATTRIBUTE_UNUSED,
>   			       const REAL_VALUE_TYPE &lh_ub ATTRIBUTE_UNUSED,
>   			       const REAL_VALUE_TYPE &rh_lb ATTRIBUTE_UNUSED,
> -			       const REAL_VALUE_TYPE &rh_ub ATTRIBUTE_UNUSED)
> -  const
> +			       const REAL_VALUE_TYPE &rh_ub ATTRIBUTE_UNUSED,
> +			       relation_kind) const
>   {
>     lb = dconstninf;
>     ub = dconstinf;
> @@ -1868,7 +1868,8 @@ class foperator_plus : public range_oper
>   		const REAL_VALUE_TYPE &lh_lb,
>   		const REAL_VALUE_TYPE &lh_ub,
>   		const REAL_VALUE_TYPE &rh_lb,
> -		const REAL_VALUE_TYPE &rh_ub) const final override
> +		const REAL_VALUE_TYPE &rh_ub,
> +		relation_kind) const final override
>     {
>       frange_arithmetic (PLUS_EXPR, type, lb, lh_lb, rh_lb, dconstninf);
>       frange_arithmetic (PLUS_EXPR, type, ub, lh_ub, rh_ub, dconstinf);
> @@ -1892,7 +1893,8 @@ class foperator_minus : public range_ope
>   		const REAL_VALUE_TYPE &lh_lb,
>   		const REAL_VALUE_TYPE &lh_ub,
>   		const REAL_VALUE_TYPE &rh_lb,
> -		const REAL_VALUE_TYPE &rh_ub) const final override
> +		const REAL_VALUE_TYPE &rh_ub,
> +		relation_kind) const final override
>     {
>       frange_arithmetic (MINUS_EXPR, type, lb, lh_lb, rh_ub, dconstninf);
>       frange_arithmetic (MINUS_EXPR, type, ub, lh_ub, rh_lb, dconstinf);
> @@ -1908,6 +1910,123 @@ class foperator_minus : public range_ope
>     }
>   } fop_minus;
>   
> +/* Wrapper around frange_arithmetics, that computes the result
> +   if inexact rounded to both directions.  Also, if one of the
> +   operands is +-0.0 and another +-INF, return +-0.0 rather than
> +   NAN.  */

s/frange_arithmetics/frange_arithmetic/

Also, would you mind written a little blurb about why it's necessary not 
to compute INF*0.0 as NAN.  I assume it's because you're using it for 
the cross product and you'll set maybe_nan separately, but it's nice to 
spell it out.

> +
> +static void
> +frange_mult (tree type, REAL_VALUE_TYPE &result_lb, REAL_VALUE_TYPE &result_ub,
> +	     const REAL_VALUE_TYPE &op1, const REAL_VALUE_TYPE &op2)
> +{
> +  if (real_iszero (&op1) && real_isinf (&op2))
> +    {
> +      result_lb = op1;
> +      if (real_isneg (&op2))
> +	result_lb = real_value_negate (&result_lb);
> +      result_ub = result_lb;
> +    }
> +  else if (real_isinf (&op1) && real_iszero (&op2))
> +    {
> +      result_lb = op2;
> +      if (real_isneg (&op1))
> +	result_lb = real_value_negate (&result_lb);
> +      result_ub = result_lb;
> +    }
> +  else
> +    {
> +      frange_arithmetic (MULT_EXPR, type, result_lb, op1, op2, dconstninf);
> +      frange_arithmetic (MULT_EXPR, type, result_ub, op1, op2, dconstinf);
> +    }
> +}
> +
> +class foperator_mult : public range_operator_float
> +{
> +  void rv_fold (REAL_VALUE_TYPE &lb, REAL_VALUE_TYPE &ub, bool &maybe_nan,
> +		tree type,
> +		const REAL_VALUE_TYPE &lh_lb,
> +		const REAL_VALUE_TYPE &lh_ub,
> +		const REAL_VALUE_TYPE &rh_lb,
> +		const REAL_VALUE_TYPE &rh_ub,
> +		relation_kind kind) const final override
> +  {
> +    REAL_VALUE_TYPE cp[8];
> +    bool is_square
> +      = (kind == VREL_EQ
> +	 && real_equal (&lh_lb, &rh_lb)
> +	 && real_equal (&lh_ub, &rh_ub)
> +	 && real_isneg (&lh_lb) == real_isneg (&rh_lb)
> +	 && real_isneg (&lh_ub) == real_isneg (&rh_ub));
> +    // Do a cross-product.
> +    frange_mult (type, cp[0], cp[4], lh_lb, rh_lb);
> +    if (is_square)
> +      {
> +	// For x * x we can just do max (lh_lb * lh_lb, lh_ub * lh_ub)
> +	// as maximum and -0.0 as minimum if 0.0 is in the range,
> +	// otherwise min (lh_lb * lh_lb, lh_ub * lh_ub).
> +	// -0.0 rather than 0.0 because VREL_EQ doesn't prove that
> +	// x and y are bitwise equal, just that they compare equal.
> +	if (real_compare (LE_EXPR, &lh_lb, &dconst0)
> +	    && real_compare (GE_EXPR, &lh_ub, &dconst0))
> +	  cp[1] = real_value_negate (&dconst0);
> +	else
> +	  cp[1] = cp[0];
> +	cp[2] = cp[0];
> +	cp[5] = cp[4];
> +	cp[6] = cp[4];
> +      }
> +    else
> +      {
> +	frange_mult (type, cp[1], cp[5], lh_lb, rh_ub);
> +	frange_mult (type, cp[2], cp[6], lh_ub, rh_lb);
> +      }
> +    frange_mult (type, cp[3], cp[7], lh_ub, rh_ub);
> +    for (int i = 1; i < 4; ++i)
> +      {
> +	if (real_less (&cp[i], &cp[0])
> +	    || (real_iszero (&cp[0]) && real_isnegzero (&cp[i])))
> +	  std::swap (cp[i], cp[0]);
> +	if (real_less (&cp[4], &cp[i + 4])
> +	    || (real_isnegzero (&cp[4]) && real_iszero (&cp[i + 4])))
> +	  std::swap (cp[i + 4], cp[4]);
> +      }
> +    lb = cp[0];
> +    ub = cp[4];
> +
> +    // If both operands are the same, then we know it can be +-0.0, or +-INF,
> +    // but not both at the same time, so it will never be invalid unless
> +    // operand was already NAN.
> +    if (is_square)
> +      maybe_nan = false;
> +    // [+-0, +-0] * [+INF,+INF] (or [-INF,-INF] or swapped is a known NAN.
> +    else if ((real_iszero (&lh_lb)
> +	      && real_iszero (&lh_ub)
> +	      && real_isinf (&rh_lb)
> +	      && real_isinf (&rh_ub, real_isneg (&rh_lb)))
> +	     || (real_iszero (&rh_lb)
> +		 && real_iszero (&rh_ub)
> +		 && real_isinf (&lh_lb)
> +		 && real_isinf (&lh_ub, real_isneg (&lh_lb))))
> +      {
> +	real_nan (&lb, "", 0, TYPE_MODE (type));
> +	ub = lb;
> +	maybe_nan = true;
> +      }
> +    // Otherwise, if one range includes zero and the other ends with +-INF,
> +    // it is a maybe NAN.
> +    else if (real_compare (LE_EXPR, &lh_lb, &dconst0)
> +	     && real_compare (GE_EXPR, &lh_ub, &dconst0)
> +	     && (real_isinf (&rh_lb) || real_isinf (&rh_ub)))
> +      maybe_nan = true;
> +    else if (real_compare (LE_EXPR, &rh_lb, &dconst0)
> +	     && real_compare (GE_EXPR, &rh_ub, &dconst0)
> +	     && (real_isinf (&lh_lb) || real_isinf (&lh_ub)))
> +      maybe_nan = true;
> +    else
> +      maybe_nan = false;
> +  }
> +} fop_mult;
> +
>   // Instantiate a range_op_table for floating point operations.
>   static floating_op_table global_floating_table;
>   
> @@ -1942,6 +2061,7 @@ floating_op_table::floating_op_table ()
>     set (NEGATE_EXPR, fop_negate);
>     set (PLUS_EXPR, fop_plus);
>     set (MINUS_EXPR, fop_minus);
> +  set (MULT_EXPR, fop_mult);
>   }
>   
>   // Return a pointer to the range_operator_float instance, if there is

It'd be nice to have some testcases.  For example, from what I can see, 
the original integer multiplication code came with some tests in 
gcc.dg/tree-ssa/vrp13.c (commit 9983270bec0a18).  It'd be nice to have 
some sanity checks, especially because so many things can go wrong with 
floats.

I'll leave it to you to decide what tests to include.

Otherwise, LGTM.

Thanks.
Aldy

On Thu, Nov 10, 2022 at 03:50:47PM +0100, Aldy Hernandez wrote:
> > @@ -1908,6 +1910,123 @@ class foperator_minus : public range_ope
> >     }
> >   } fop_minus;
> > +/* Wrapper around frange_arithmetics, that computes the result
> > +   if inexact rounded to both directions.  Also, if one of the
> > +   operands is +-0.0 and another +-INF, return +-0.0 rather than
> > +   NAN.  */
> 
> s/frange_arithmetics/frange_arithmetic/
> 
> Also, would you mind written a little blurb about why it's necessary not to
> compute INF*0.0 as NAN.  I assume it's because you're using it for the cross
> product and you'll set maybe_nan separately, but it's nice to spell it out.

This made me think about it some more and I'll need to play around with it
some more, perhaps the right thing is similarly to what I've attached for
division to handle special cases upfront and call frange_arithmetic only
for the safe cases.
E.g. one case which the posted foperator_mult handles pessimistically is
[0.0, 10.0] * [INF, INF].  This should be just [INF, INF] +-NAN IMHO,
because the 0.0 * INF case will result in NAN, while
nextafter (0.0, 1.0) * INF
will be already INF and everything larger as well.
I could in frange_mult be very conservative and for the 0 * INF cases
set result_lb and result_ub to [0.0, INF] range (corresponding signs
depending on the xor of sign of ops), but that would be quite pessimistic as
well.  If one has:
[0.0, 0.0] * [10.0, INF], the result should be just [0.0, 0.0] +-NAN,
because again 0.0 * INF is NAN, but 0.0 * nextafter (INF, 0.0) is already 0.0.

Note, the is_square case doesn't suffer from all of this mess, the result
is never NAN (unless operand is NAN).

> It'd be nice to have some testcases.  For example, from what I can see, the
> original integer multiplication code came with some tests in
> gcc.dg/tree-ssa/vrp13.c (commit 9983270bec0a18).  It'd be nice to have some
> sanity checks, especially because so many things can go wrong with floats.
> 
> I'll leave it to you to decide what tests to include.

I've tried following, but it suffers from various issues:
1) we don't handle __builtin_signbit (whatever) == 0 (or != 0) as guarantee
   that in the guarded code whatever has signbit 0 or 1
2) __builtin_isinf (x) > 0 is lowered to x > DBL_MAX, but unfortunately we don't
   infer from that [INF,INF] range, but [DBL_MAX, INF] range
3) what I wrote above, I think we don't handle [0, 2] * [INF, INF] right but
   due to 2) we can't see it

So, maybe for now a selftest will be better than a testcase, or
alternatively a plugin test which acts like a selftest.

/* { dg-do compile { target { ! { vax-*-* powerpc-*-*spe pdp11-*-* } } } } */
/* { dg-options "-O2 -fno-trapping-math -fno-signaling-nans -fsigned-zeros -fno-tree-fre -fno-tree-dominator-opts -fno-thread-jumps -fdump-tree-optimized" } */
/* { dg-add-options ieee } */

void
foo (double x, double y)
{
  const double inf = __builtin_inf ();
  const double minf = -inf;
  if (__builtin_isnan (x) || __builtin_isnan (y))
    return;
#define TEST(n, xl, xu, yl, yu, rl, ru, nan) \
  if ((__builtin_isinf (xl) > 0				\
       ? x > 0.0 && __builtin_isinf (x)			\
       : __builtin_isinf (xu) < 0			\
       ? x < 0.0 && __builtin_isinf (x)			\
       : x >= xl && x <= xu				\
	 && (xl != 0.0					\
	     || __builtin_signbit (xl)			\
	     || !__builtin_signbit (x))			\
	 && (xu != 0.0					\
	     || !__builtin_signbit (xu)			\
	     || __builtin_signbit (x)))			\
      && (__builtin_isinf (yl) > 0			\
	  ? y > 0.0 && __builtin_isinf (y)		\
	  : __builtin_isinf (yu) < 0			\
	  ? y < 0.0 && __builtin_isinf (y)		\
	  : y >= yl && y <= yu				\
	    && (yl != 0.0				\
		|| __builtin_signbit (yl)		\
		|| !__builtin_signbit (y))		\
	    && (yu != 0.0				\
		|| !__builtin_signbit (yu)		\
		|| __builtin_signbit (y))))		\
    {							\
      double r##n = x * y;				\
      if (nan == 2)					\
	{						\
	  if (!__builtin_isnan (r##n))			\
	    __builtin_abort ();				\
	}						\
      else if (nan == 1)				\
	{						\
	  if (!__builtin_isnan (r##n))			\
	    {						\
	      if (r##n < rl || r##n > ru)		\
		__builtin_abort ();			\
	    }						\
	}						\
      else						\
	{						\
	  if (__builtin_isnan (r##n))			\
	    __builtin_abort ();				\
	  if (r##n < rl || r##n > ru)			\
	    __builtin_abort ();				\
	}						\
    }
#define TEST2(n, xl, xu, rl, ru) \
  if (__builtin_isinf (xl) > 0				\
      ? x > 0.0 && __builtin_isinf (x)			\
      : __builtin_isinf (xu) < 0			\
      ? x < 0.0 && __builtin_isinf (x)			\
      : x >= xl && x <= xu				\
	&& (xl != 0.0					\
	    || __builtin_signbit (xl)			\
	    || !__builtin_signbit (x))			\
	&& (xu != 0.0					\
	    || !__builtin_signbit (xu)			\
	    || __builtin_signbit (x)))			\
    {							\
      double s##n = x * x;				\
      if (__builtin_isnan (s##n))			\
	__builtin_abort ();				\
      if (s##n < rl || s##n > ru)			\
	__builtin_abort ();				\
    }
  TEST (1,	2.0, 4.0,	6.0, 8.0,	12.0, 32.0, 0);
  TEST (2,	-2.0, 3.0,	-7.0, 4.0,	-21.0, 14.0, 0);
  TEST (3,	-9.0, 5.0,	8.0, 10.0,	-90.0, 50.0, 0);
  TEST (4,	-0.0, 0.0,	16.0, 32.0,	-0.0, 0.0, 0);
  TEST (5,	-0.0, 0.0,	16.0, 32.0,	-0.0, 0.0, 0);
  TEST (6,	0.0, 0.0,	16.0, 32.0,	0.0, 0.0, 0);
  TEST (7,	0.0, 0.0,	16.0, 32.0,	0.0, 0.0, 0);
  TEST (8,	-0.0, -0.0,	16.0, 32.0,	-0.0, -0.0, 0);
  TEST (9,	-0.0, -0.0,	16.0, 32.0,	-0.0, -0.0, 0);
  TEST (10,	minf, inf,	minf, inf,	minf, inf, 1);
  TEST (11,	-0.0, 0.0,	128.0, inf,	-0.0, 0.0, 1);
  TEST (12,	-0.0, 0.0,	inf, inf,	0.0, 0.0, 2);
  TEST (13,	minf, minf,	0.0, 0.0,	0.0, 0.0, 2);
  TEST (14,	0.0, 2.0,	inf, inf,	inf, inf, 1);
  TEST (15,	0.0, 2.0,	minf, minf,	minf, minf, 1);
  TEST (16,	inf, inf,	-0.0, 2.0,	minf, inf, 1);
  TEST (17,	minf, minf,	-0.0, 3.0,	minf, inf, 1);

  TEST2 (1,	2.0, 4.0,			4.0, 16.0);
  TEST2 (2,	-0.0, 0.0,			-0.0, 0.0);
  TEST2 (3,	0.0, inf,			0.0, inf);
  TEST2 (4,	inf, inf,			inf, inf);
}


	Jakub

On Thu, Nov 10, 2022 at 08:20:06PM +0100, Jakub Jelinek via Gcc-patches wrote:
> So, maybe for now a selftest will be better than a testcase, or
> alternatively a plugin test which acts like a selftest.

For a testsuite/g*.dg/plugin/ range-op-plugin.c test, would be
nice to write a short plugin that does:
1) register 2 new attributes, say gcc_range and gcc_expected_range,
   parse their arguments
2) registers some new pass (dunno where it would be best, say before evrp
   or so), which would:
   - for function parameters with gcc_range attributes set global? range
     on default def of that parameter
   - if function itself has a gcc_expected_range attribute, propagate
     ranges from arguments to the function result and compare against
     gcc_expected_range
   - if function itself has gcc_range attribute and one of the arguments
     gcc_expected_range itself, try to propagate range backwards from
     result to the argument
Then we could say write tests like:
__attribute__((gcc_expected_range (12.0, 32.0, 0))) double
foo (double x __attribute__((gcc_range (2.0, 4.0, 0))), double y __attribute__((gcc_range (6.0, 8.0, 0))))
{
  return x * y;
}
with for floating point types (parameter type or function result type)
the arguments of the attribute being (constant folded)
low bound, high bound, integer about NAN (say 0 meaning clear_nan,
bit 0 meaning +NAN, bit 1 -NAN, bit 2 meaning known NAN (then
the 2 bounds would be ignored)).
Eventually we could do something similar for integral types, pointer types
etc.
I think this would be far more useful compared to writing selftests for it,
and compared to the testcase I've posted would be easier to request or check
exact range rather than a rough range.  And we could easily
test not just very simple binary ops (in both directions), but builtin calls
etc.

Thoughts on this?

I can try to write a plugin that registers the attributes, parses them
and registers a pass, but would appreciate your or Andrew's help in filling
the actual pass.

	Jakub

On 11/10/22 20:20, Jakub Jelinek wrote:
> On Thu, Nov 10, 2022 at 03:50:47PM +0100, Aldy Hernandez wrote:
>>> @@ -1908,6 +1910,123 @@ class foperator_minus : public range_ope
>>>      }
>>>    } fop_minus;
>>> +/* Wrapper around frange_arithmetics, that computes the result
>>> +   if inexact rounded to both directions.  Also, if one of the
>>> +   operands is +-0.0 and another +-INF, return +-0.0 rather than
>>> +   NAN.  */
>>
>> s/frange_arithmetics/frange_arithmetic/
>>
>> Also, would you mind written a little blurb about why it's necessary not to
>> compute INF*0.0 as NAN.  I assume it's because you're using it for the cross
>> product and you'll set maybe_nan separately, but it's nice to spell it out.
> 
> This made me think about it some more and I'll need to play around with it
> some more, perhaps the right thing is similarly to what I've attached for
> division to handle special cases upfront and call frange_arithmetic only
> for the safe cases.
> E.g. one case which the posted foperator_mult handles pessimistically is
> [0.0, 10.0] * [INF, INF].  This should be just [INF, INF] +-NAN IMHO,
> because the 0.0 * INF case will result in NAN, while
> nextafter (0.0, 1.0) * INF
> will be already INF and everything larger as well.
> I could in frange_mult be very conservative and for the 0 * INF cases
> set result_lb and result_ub to [0.0, INF] range (corresponding signs
> depending on the xor of sign of ops), but that would be quite pessimistic as
> well.  If one has:
> [0.0, 0.0] * [10.0, INF], the result should be just [0.0, 0.0] +-NAN,
> because again 0.0 * INF is NAN, but 0.0 * nextafter (INF, 0.0) is already 0.0.
> 
> Note, the is_square case doesn't suffer from all of this mess, the result
> is never NAN (unless operand is NAN).
> 
>> It'd be nice to have some testcases.  For example, from what I can see, the
>> original integer multiplication code came with some tests in
>> gcc.dg/tree-ssa/vrp13.c (commit 9983270bec0a18).  It'd be nice to have some
>> sanity checks, especially because so many things can go wrong with floats.
>>
>> I'll leave it to you to decide what tests to include.
> 
> I've tried following, but it suffers from various issues:
> 1) we don't handle __builtin_signbit (whatever) == 0 (or != 0) as guarantee
>     that in the guarded code whatever has signbit 0 or 1

We have a range-op entry for __builtin_signbit in cfn_signbit.  Is this 
a shortcoming of this code, or something else?

> 2) __builtin_isinf (x) > 0 is lowered to x > DBL_MAX, but unfortunately we don't
>     infer from that [INF,INF] range, but [DBL_MAX, INF] range
> 3) what I wrote above, I think we don't handle [0, 2] * [INF, INF] right but
>     due to 2) we can't see it

Doesn't this boil down to a representation issue?  I wonder if we should 
bite the bullet and tweak build_gt() and build_lt() to represent open 
ranges.  In theory it should be one more/less ULP, while adjusting for 
HONOR_INFINITIES.

If the signbit issue were resolved and we could represent > and < 
properly, would that allow us to write proper testcases without having 
to writing a plug-in (which I assume is a lot harder)?

Aldy

> 
> So, maybe for now a selftest will be better than a testcase, or
> alternatively a plugin test which acts like a selftest.
> 
> /* { dg-do compile { target { ! { vax-*-* powerpc-*-*spe pdp11-*-* } } } } */
> /* { dg-options "-O2 -fno-trapping-math -fno-signaling-nans -fsigned-zeros -fno-tree-fre -fno-tree-dominator-opts -fno-thread-jumps -fdump-tree-optimized" } */
> /* { dg-add-options ieee } */
> 
> void
> foo (double x, double y)
> {
>    const double inf = __builtin_inf ();
>    const double minf = -inf;
>    if (__builtin_isnan (x) || __builtin_isnan (y))
>      return;
> #define TEST(n, xl, xu, yl, yu, rl, ru, nan) \
>    if ((__builtin_isinf (xl) > 0				\
>         ? x > 0.0 && __builtin_isinf (x)			\
>         : __builtin_isinf (xu) < 0			\
>         ? x < 0.0 && __builtin_isinf (x)			\
>         : x >= xl && x <= xu				\
> 	 && (xl != 0.0					\
> 	     || __builtin_signbit (xl)			\
> 	     || !__builtin_signbit (x))			\
> 	 && (xu != 0.0					\
> 	     || !__builtin_signbit (xu)			\
> 	     || __builtin_signbit (x)))			\
>        && (__builtin_isinf (yl) > 0			\
> 	  ? y > 0.0 && __builtin_isinf (y)		\
> 	  : __builtin_isinf (yu) < 0			\
> 	  ? y < 0.0 && __builtin_isinf (y)		\
> 	  : y >= yl && y <= yu				\
> 	    && (yl != 0.0				\
> 		|| __builtin_signbit (yl)		\
> 		|| !__builtin_signbit (y))		\
> 	    && (yu != 0.0				\
> 		|| !__builtin_signbit (yu)		\
> 		|| __builtin_signbit (y))))		\
>      {							\
>        double r##n = x * y;				\
>        if (nan == 2)					\
> 	{						\
> 	  if (!__builtin_isnan (r##n))			\
> 	    __builtin_abort ();				\
> 	}						\
>        else if (nan == 1)				\
> 	{						\
> 	  if (!__builtin_isnan (r##n))			\
> 	    {						\
> 	      if (r##n < rl || r##n > ru)		\
> 		__builtin_abort ();			\
> 	    }						\
> 	}						\
>        else						\
> 	{						\
> 	  if (__builtin_isnan (r##n))			\
> 	    __builtin_abort ();				\
> 	  if (r##n < rl || r##n > ru)			\
> 	    __builtin_abort ();				\
> 	}						\
>      }
> #define TEST2(n, xl, xu, rl, ru) \
>    if (__builtin_isinf (xl) > 0				\
>        ? x > 0.0 && __builtin_isinf (x)			\
>        : __builtin_isinf (xu) < 0			\
>        ? x < 0.0 && __builtin_isinf (x)			\
>        : x >= xl && x <= xu				\
> 	&& (xl != 0.0					\
> 	    || __builtin_signbit (xl)			\
> 	    || !__builtin_signbit (x))			\
> 	&& (xu != 0.0					\
> 	    || !__builtin_signbit (xu)			\
> 	    || __builtin_signbit (x)))			\
>      {							\
>        double s##n = x * x;				\
>        if (__builtin_isnan (s##n))			\
> 	__builtin_abort ();				\
>        if (s##n < rl || s##n > ru)			\
> 	__builtin_abort ();				\
>      }
>    TEST (1,	2.0, 4.0,	6.0, 8.0,	12.0, 32.0, 0);
>    TEST (2,	-2.0, 3.0,	-7.0, 4.0,	-21.0, 14.0, 0);
>    TEST (3,	-9.0, 5.0,	8.0, 10.0,	-90.0, 50.0, 0);
>    TEST (4,	-0.0, 0.0,	16.0, 32.0,	-0.0, 0.0, 0);
>    TEST (5,	-0.0, 0.0,	16.0, 32.0,	-0.0, 0.0, 0);
>    TEST (6,	0.0, 0.0,	16.0, 32.0,	0.0, 0.0, 0);
>    TEST (7,	0.0, 0.0,	16.0, 32.0,	0.0, 0.0, 0);
>    TEST (8,	-0.0, -0.0,	16.0, 32.0,	-0.0, -0.0, 0);
>    TEST (9,	-0.0, -0.0,	16.0, 32.0,	-0.0, -0.0, 0);
>    TEST (10,	minf, inf,	minf, inf,	minf, inf, 1);
>    TEST (11,	-0.0, 0.0,	128.0, inf,	-0.0, 0.0, 1);
>    TEST (12,	-0.0, 0.0,	inf, inf,	0.0, 0.0, 2);
>    TEST (13,	minf, minf,	0.0, 0.0,	0.0, 0.0, 2);
>    TEST (14,	0.0, 2.0,	inf, inf,	inf, inf, 1);
>    TEST (15,	0.0, 2.0,	minf, minf,	minf, minf, 1);
>    TEST (16,	inf, inf,	-0.0, 2.0,	minf, inf, 1);
>    TEST (17,	minf, minf,	-0.0, 3.0,	minf, inf, 1);
> 
>    TEST2 (1,	2.0, 4.0,			4.0, 16.0);
>    TEST2 (2,	-0.0, 0.0,			-0.0, 0.0);
>    TEST2 (3,	0.0, inf,			0.0, inf);
>    TEST2 (4,	inf, inf,			inf, inf);
> }
> 
> 
> 	Jakub
>

On Fri, Nov 11, 2022 at 11:01:38AM +0100, Aldy Hernandez wrote:
> > I've tried following, but it suffers from various issues:
> > 1) we don't handle __builtin_signbit (whatever) == 0 (or != 0) as guarantee
> >     that in the guarded code whatever has signbit 0 or 1
> 
> We have a range-op entry for __builtin_signbit in cfn_signbit.  Is this a
> shortcoming of this code, or something else?

Dunno, I admit I haven't investigated it much.  I just saw it when putting
a breakpoint on the mult fold_range.

> > 2) __builtin_isinf (x) > 0 is lowered to x > DBL_MAX, but unfortunately we don't
> >     infer from that [INF,INF] range, but [DBL_MAX, INF] range
> > 3) what I wrote above, I think we don't handle [0, 2] * [INF, INF] right but
> >     due to 2) we can't see it
> 
> Doesn't this boil down to a representation issue?  I wonder if we should
> bite the bullet and tweak build_gt() and build_lt() to represent open
> ranges.  In theory it should be one more/less ULP, while adjusting for
> HONOR_INFINITIES.

At least with the exception of MODE_COMPOSITE_P, I think we don't need
to introduce open ranges (and who cares about MODE_COMPOSITE_P if it is
conservatively correct).
For other floats, I think
x > cst
is always equivalent to
x >= nextafter (cst, inf)
and
x < cst
is always equivalent to
x <= nextafter (cst, -inf)
except for the signed zero cases which needs tiny bit more thought.
So, if we have
if (x > DBL_MAX)
then in code guarded by that we can just use [INF, INF] as range.

> If the signbit issue were resolved and we could represent > and < properly,
> would that allow us to write proper testcases without having to writing a
> plug-in (which I assume is a lot harder)?

I don't know, we'd need to see.
First work out on all the issues that result on the testcase the operand
ranges aren't exactly what we want (whether on the testcase side or on the
range-ops side or wherever) and once that looks ok, see if the ranges
on the rN/sN vars are correct and if so, watch what hasn't been folded away
and why.
I think the plugin would be 100-200 lines of code and then we could just
write multiple testcases against the plugin.

	Jakub

On 11/11/22 11:47, Jakub Jelinek wrote:
> On Fri, Nov 11, 2022 at 11:01:38AM +0100, Aldy Hernandez wrote:
>>> I've tried following, but it suffers from various issues:
>>> 1) we don't handle __builtin_signbit (whatever) == 0 (or != 0) as guarantee
>>>      that in the guarded code whatever has signbit 0 or 1
>>
>> We have a range-op entry for __builtin_signbit in cfn_signbit.  Is this a
>> shortcoming of this code, or something else?
> 
> Dunno, I admit I haven't investigated it much.  I just saw it when putting
> a breakpoint on the mult fold_range.

Could you send me a small testcase.  I can look into that.

> 
>>> 2) __builtin_isinf (x) > 0 is lowered to x > DBL_MAX, but unfortunately we don't
>>>      infer from that [INF,INF] range, but [DBL_MAX, INF] range
>>> 3) what I wrote above, I think we don't handle [0, 2] * [INF, INF] right but
>>>      due to 2) we can't see it
>>
>> Doesn't this boil down to a representation issue?  I wonder if we should
>> bite the bullet and tweak build_gt() and build_lt() to represent open
>> ranges.  In theory it should be one more/less ULP, while adjusting for
>> HONOR_INFINITIES.
> 
> At least with the exception of MODE_COMPOSITE_P, I think we don't need
> to introduce open ranges (and who cares about MODE_COMPOSITE_P if it is
> conservatively correct).
> For other floats, I think
> x > cst
> is always equivalent to
> x >= nextafter (cst, inf)
> and
> x < cst
> is always equivalent to
> x <= nextafter (cst, -inf)
> except for the signed zero cases which needs tiny bit more thought.
> So, if we have
> if (x > DBL_MAX)
> then in code guarded by that we can just use [INF, INF] as range.

Yeah, yeah.  That's exactly what I meant... using nextafter.  I'll look 
into that, as there seems there's more than one issue related to our 
lack of precision in representing < and >.

> 
>> If the signbit issue were resolved and we could represent > and < properly,
>> would that allow us to write proper testcases without having to writing a
>> plug-in (which I assume is a lot harder)?
> 
> I don't know, we'd need to see.
> First work out on all the issues that result on the testcase the operand
> ranges aren't exactly what we want (whether on the testcase side or on the
> range-ops side or wherever) and once that looks ok, see if the ranges
> on the rN/sN vars are correct and if so, watch what hasn't been folded away
> and why.
> I think the plugin would be 100-200 lines of code and then we could just
> write multiple testcases against the plugin.

If you think the plug-in will get better test coverage, by all means.  I 
was just trying to save you/us some work.

Andrew, do you have any thoughts on the plug-in?

Aldy