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

[Jakub Jelinek <jakub@redhat.com>]

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