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

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

[Jakub Jelinek]

Hi!

Here is the floating point division fold_range implementation,
as I wrote in the last mail, we could outline some of the common parts
into static methods with descriptive names and share them between
foperator_div and foperator_mult.

Bootstrapped/regtested on top of the earlier version of the multiplication
fold_range on x86_64-linux and i686-linux, regressions are
+FAIL: gcc.dg/pr95115.c execution test
+FAIL: libphobos.phobos/std/math/hardware.d execution test
+FAIL: libphobos.phobos_shared/std/math/hardware.d execution test
The first test is we have:
  # RANGE [frange] double [] +-NAN
  _3 =  Inf /  Inf;
  if (_3 ord _3)
    goto <bb 3>; [INV]
  else
    goto <bb 4>; [INV]

  <bb 3> :
  abort ();

  <bb 4> :
before evrp, the range is correct, Inf / Inf is known NAN of unknown
sign.  evrp correctly folds _3 ord _3 into false and the
  _3 =  Inf /  Inf;
remains in the IL, but then comes dse1 and removes it as dead
statement.  So, I think yet another example of the PR107608 problems
where DCE? removes dead statements which raise floating point exceptions.
And -fno-delete-dead-exceptions doesn't help.

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

	PR tree-optimization/107569
	* range-op-float.cc (foperator_div): New class.
	(floating_op_table::floating_op_table): Use foperator_div
	for RDIV_EXPR.

--- gcc/range-op-float.cc.jj	2022-11-10 12:31:57.987917289 +0100
+++ gcc/range-op-float.cc	2022-11-10 17:04:35.743056880 +0100
@@ -2027,6 +2027,183 @@ class foperator_mult : public range_oper
   }
 } fop_mult;
 
+class foperator_div : 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) const final override
+  {
+    // +-0.0 / +-0.0 or +-INF / +-INF is a known NAN.
+    if ((real_iszero (&lh_lb)
+	 && real_iszero (&lh_ub)
+	 && real_iszero (&rh_lb)
+	 && real_iszero (&rh_ub))
+	|| (real_isinf (&lh_lb)
+	    && real_isinf (&lh_ub, real_isneg (&lh_lb))
+	    && real_isinf (&rh_lb)
+	    && real_isinf (&rh_ub, real_isneg (&rh_lb))))
+      {
+	real_nan (&lb, "", 0, TYPE_MODE (type));
+	ub = lb;
+	maybe_nan = true;
+	return;
+      }
+
+    bool both_maybe_zero = false;
+    bool both_maybe_inf = false;
+    bool must_have_signbit_zero = false;
+    bool must_have_signbit_nonzero = false;
+
+    // If +-0.0 is in both ranges, it is a maybe NAN.
+    if (real_compare (LE_EXPR, &lh_lb, &dconst0)
+	&& real_compare (GE_EXPR, &lh_ub, &dconst0)
+	&& real_compare (LE_EXPR, &rh_lb, &dconst0)
+	&& real_compare (GE_EXPR, &rh_ub, &dconst0))
+      {
+	both_maybe_zero = true;
+	maybe_nan = true;
+      }
+    // If +-INF is in both ranges, it is a maybe NAN.
+    else if ((real_isinf (&lh_lb) || real_isinf (&lh_ub))
+	     && (real_isinf (&rh_lb) || real_isinf (&rh_ub)))
+      {
+	both_maybe_inf = true;
+	maybe_nan = true;
+      }
+    else
+      maybe_nan = false;
+
+    if (real_isneg (&lh_lb) == real_isneg (&lh_ub)
+	&& real_isneg (&rh_lb) == real_isneg (&rh_ub))
+      {
+	if (real_isneg (&lh_lb) == real_isneg (&rh_ub))
+	  must_have_signbit_zero = true;
+	else
+	  must_have_signbit_nonzero = true;
+      }
+
+    // If dividend must be zero, the range is just +-0
+    // (including if the divisor is +-INF).
+    // If divisor must be +-INF, the range is just +-0
+    // (including if the dividend is zero).
+    if ((real_iszero (&lh_lb) && real_iszero (&lh_ub))
+	|| real_isinf (&rh_lb, false)
+	|| real_isinf (&rh_ub, true))
+      {
+	ub = lb = dconst0;
+	// If all the boundary signs are the same, [+0.0, +0.0].
+	if (must_have_signbit_zero)
+	  ;
+	// If divisor and dividend must have different signs,
+	// [-0.0, -0.0].
+	else if (must_have_signbit_nonzero)
+	  ub = lb = real_value_negate (&dconst0);
+	// Otherwise -> [-0.0, +0.0].
+	else
+	  lb = real_value_negate (&dconst0);
+	return;
+      }
+
+    // If divisor must be zero, the range is just +-INF
+    // (including if the dividend is +-INF).
+    // If dividend must be +-INF, the range is just +-INF
+    // (including if the dividend is zero).
+    if ((real_iszero (&rh_lb) && real_iszero (&rh_ub))
+	|| real_isinf (&lh_lb, false)
+	|| real_isinf (&lh_ub, true))
+      {
+	// If all the boundary signs are the same, [+INF, +INF].
+	if (must_have_signbit_zero)
+	  ub = lb = dconstinf;
+	// If divisor and dividend must have different signs,
+	// [-INF, -INF].
+	else if (must_have_signbit_nonzero)
+	  ub = lb = dconstninf;
+	// Otherwise -> [-INF, +INF] (-INF or +INF).
+	else
+	  {
+	    lb = dconstninf;
+	    ub = dconstinf;
+	  }
+	return;
+      }
+
+    // Otherwise if both operands may be zero, divisor could be
+    // nextafter(0.0, +-1.0) and dividend +-0.0
+    // in which case result is going to INF or vice versa and
+    // result +0.0.  So, all we can say for that case is if the
+    // signs of divisor and dividend are always the same we have
+    // [+0.0, +INF], if they are always different we have
+    // [-INF, -0.0].  If they vary, VARING.
+    // If both may be +-INF, divisor could be INF and dividend FLT_MAX,
+    // in which case result is going to INF or vice versa and
+    // result +0.0.  So, all we can say for that case is if the
+    // signs of divisor and dividend are always the same we have
+    // [+0.0, +INF], if they are always different we have
+    // [-INF, -0.0].  If they vary, VARYING.
+    if (both_maybe_zero || both_maybe_inf)
+      {
+	if (must_have_signbit_zero)
+	  {
+	    lb = dconst0;
+	    ub = dconstinf;
+	  }
+	else if (must_have_signbit_nonzero)
+	  {
+	    lb = dconstninf;
+	    ub = real_value_negate (&dconst0);
+	  }
+	else
+	  {
+	    lb = dconstninf;
+	    ub = dconstinf;
+	  }
+	return;
+      }
+
+    REAL_VALUE_TYPE cp[8];
+    // Do a cross-division.  At this point none of the divisions should
+    // produce a NAN.
+    gcc_assert (!maybe_nan);
+    frange_arithmetic (RDIV_EXPR, type, cp[0], lh_lb, rh_lb, dconstninf);
+    frange_arithmetic (RDIV_EXPR, type, cp[1], lh_lb, rh_ub, dconstninf);
+    frange_arithmetic (RDIV_EXPR, type, cp[2], lh_ub, rh_lb, dconstninf);
+    frange_arithmetic (RDIV_EXPR, type, cp[3], lh_ub, rh_ub, dconstninf);
+    frange_arithmetic (RDIV_EXPR, type, cp[4], lh_lb, rh_lb, dconstinf);
+    frange_arithmetic (RDIV_EXPR, type, cp[5], lh_lb, rh_ub, dconstinf);
+    frange_arithmetic (RDIV_EXPR, type, cp[6], lh_ub, rh_lb, dconstinf);
+    frange_arithmetic (RDIV_EXPR, type, cp[7], lh_ub, rh_ub, dconstinf);
+
+    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 divisor may be zero (but is not known to be only zero),
+    // and dividend can't be zero, the range can go up to -INF or +INF
+    // depending on the signs.
+    if (real_compare (LE_EXPR, &rh_lb, &dconst0)
+	&& real_compare (GE_EXPR, &rh_ub, &dconst0))
+      {
+	if (!must_have_signbit_zero)
+	  real_inf (&lb, true);
+	if (!must_have_signbit_nonzero)
+	  real_inf (&ub, false);
+      }
+  }
+} fop_div;
+
 // Instantiate a range_op_table for floating point operations.
 static floating_op_table global_floating_table;
 
@@ -2062,6 +2239,7 @@ floating_op_table::floating_op_table ()
   set (PLUS_EXPR, fop_plus);
   set (MINUS_EXPR, fop_minus);
   set (MULT_EXPR, fop_mult);
+  set (RDIV_EXPR, fop_div);
 }
 
 // Return a pointer to the range_operator_float instance, if there is

	Jakub

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

[Jan-Benedict Glaw]

[-- Attachment #1: Type: text/plain, Size: 822 bytes --]

Hi Jakub,

On Fri, 2022-11-11 10:09:42 +0100, Jakub Jelinek via Gcc-patches <gcc-patches@gcc.gnu.org> wrote:
> Here is the floating point division fold_range implementation,
> as I wrote in the last mail, we could outline some of the common parts
> into static methods with descriptive names and share them between
> foperator_div and foperator_mult.
[...]

I'm running a slightly hacked [glibc]/scripts/build-many-glibcs.py to
to CI builds for glibc as well by now (hacked to allow for GCC master
being used) and this GCC commit
(2d5c4a16dd833aa083f13dd3e78e3ef38afe6ebb) triggers glibc's
elf/check-localplt testcase to fail, though just for
sparc64-linux-gnu. (As I just started with glibc checks, it took me a
while to realize this was a real regression and not a flaw in my
setup.)

MfG, JBG

-- 

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 195 bytes --]

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

[Joseph Myers]

On Tue, 22 Nov 2022, Jan-Benedict Glaw wrote:

> I'm running a slightly hacked [glibc]/scripts/build-many-glibcs.py to
> to CI builds for glibc as well by now (hacked to allow for GCC master
> being used) and this GCC commit
> (2d5c4a16dd833aa083f13dd3e78e3ef38afe6ebb) triggers glibc's
> elf/check-localplt testcase to fail, though just for
> sparc64-linux-gnu. (As I just started with glibc checks, it took me a
> while to realize this was a real regression and not a flaw in my
> setup.)

I think the appropriate fix is to update the relevant localplt.data (to 
add the relevant libgcc symbol marked with "?" as optional), I don't think 
there's a GCC bug here.

-- 
Joseph S. Myers
joseph@codesourcery.com

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

[Iain Buclaw]

Excerpts from Jakub Jelinek via Gcc-patches's message of November 11, 2022 10:09 am:
> Hi!
> 
> Here is the floating point division fold_range implementation,
> as I wrote in the last mail, we could outline some of the common parts
> into static methods with descriptive names and share them between
> foperator_div and foperator_mult.
> 
> Bootstrapped/regtested on top of the earlier version of the multiplication
> fold_range on x86_64-linux and i686-linux, regressions are
> +FAIL: gcc.dg/pr95115.c execution test
> +FAIL: libphobos.phobos/std/math/hardware.d execution test
> +FAIL: libphobos.phobos_shared/std/math/hardware.d execution test

I've had some time to look at the Phobos failures, and seems to me that
it's a poorly written test.

    pragma(inline, false) static void blockopt(ref real x) {}
    real a = 3.5;
    // Set all the flags to zero
    resetIeeeFlags();
    assert(!ieeeFlags.divByZero);
    blockopt(a); // avoid constant propagation by the optimizer
    // Perform a division by zero.
    a /= 0.0L;
    assert(a == real.infinity);
    assert(ieeeFlags.divByZero);
    blockopt(a); // avoid constant propagation by the optimizer


1. Since this patch, that `a /= 0.0L` operation no longer appears in the
final assembly - so no divide-by-zero flags are raised.

2. Whoever introduced blockopt() perhaps did not understand that
`a /= 0.0L` is not safe from constant propagation just because it is
surrounded by some uninlinable call.

I'll fix the test in upstream, it should really be something like:

    pragma(inline, false)
    static real forceDiv(real x, real y) { return x / y; }
    a = forceDiv(a, 0.0L);
    assert(a == real.infinity);
    assert(ieeeFlags.divByZero);


Regards,
Iain.