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

[Aldy Hernandez <aldyh@redhat.com>]

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
>