From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: from us-smtp-delivery-124.mimecast.com (us-smtp-delivery-124.mimecast.com [170.10.133.124]) by sourceware.org (Postfix) with ESMTPS id F02193858D39 for ; Fri, 11 Nov 2022 08:52:59 +0000 (GMT) DMARC-Filter: OpenDMARC Filter v1.4.1 sourceware.org F02193858D39 Authentication-Results: sourceware.org; dmarc=pass (p=none dis=none) header.from=redhat.com Authentication-Results: sourceware.org; spf=pass smtp.mailfrom=redhat.com DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=redhat.com; s=mimecast20190719; t=1668156779; h=from:from:reply-to:reply-to:subject:subject:date:date: message-id:message-id:to:to:cc:cc:mime-version:mime-version: content-type:content-type:in-reply-to:in-reply-to: references:references; bh=El5tm9M3gC3oV9T5JbZ2AgloqS1+tMhuqhGa8e3xJFY=; b=DPtI2WMn5FtcdgEgsl+xNuhwNx8iRXcl+Gu428p6dN+kpLtqnMXw1wOO0c7B5rzFTPivBe isDznvDYmxexDAEuFQ2QJaqc+EnXj2hJVbzCmSCsK6KBqLIlY59b423U2EZ5N4n1SmOGuU X0jWx1eQ638C3/MavAWr0CDUF7yzt8o= Received: from mimecast-mx02.redhat.com (mimecast-mx02.redhat.com [66.187.233.88]) by relay.mimecast.com with ESMTP with STARTTLS (version=TLSv1.2, cipher=TLS_ECDHE_RSA_WITH_AES_256_GCM_SHA384) id us-mta-466-nEwyW_xoObW0L2Wr7qKgLQ-1; Fri, 11 Nov 2022 03:52:57 -0500 X-MC-Unique: nEwyW_xoObW0L2Wr7qKgLQ-1 Received: from smtp.corp.redhat.com (int-mx06.intmail.prod.int.rdu2.redhat.com [10.11.54.6]) (using TLSv1.2 with cipher AECDH-AES256-SHA (256/256 bits)) (No client certificate requested) by mimecast-mx02.redhat.com (Postfix) with ESMTPS id 90070185A7A3 for ; Fri, 11 Nov 2022 08:52:57 +0000 (UTC) Received: from tucnak.zalov.cz (unknown [10.39.192.38]) by smtp.corp.redhat.com (Postfix) with ESMTPS id 2A9862166B36; Fri, 11 Nov 2022 08:52:57 +0000 (UTC) Received: from tucnak.zalov.cz (localhost [127.0.0.1]) by tucnak.zalov.cz (8.17.1/8.17.1) with ESMTPS id 2AB8qs4j3248960 (version=TLSv1.3 cipher=TLS_AES_256_GCM_SHA384 bits=256 verify=NOT); Fri, 11 Nov 2022 09:52:55 +0100 Received: (from jakub@localhost) by tucnak.zalov.cz (8.17.1/8.17.1/Submit) id 2AB8qspn3248959; Fri, 11 Nov 2022 09:52:54 +0100 Date: Fri, 11 Nov 2022 09:52:53 +0100 From: Jakub Jelinek To: Aldy Hernandez Cc: gcc-patches@gcc.gnu.org Subject: [PATCH] range-op, v2: Implement floating point multiplication fold_range [PR107569] Message-ID: Reply-To: Jakub Jelinek References: <7304acea-b8bb-c930-546a-db6e3b3ff96b@redhat.com> MIME-Version: 1.0 In-Reply-To: X-Scanned-By: MIMEDefang 3.1 on 10.11.54.6 X-Mimecast-Spam-Score: 0 X-Mimecast-Originator: redhat.com Content-Type: text/plain; charset=us-ascii Content-Disposition: inline X-Spam-Status: No, score=-3.5 required=5.0 tests=BAYES_00,DKIMWL_WL_HIGH,DKIM_SIGNED,DKIM_VALID,DKIM_VALID_AU,DKIM_VALID_EF,RCVD_IN_DNSWL_NONE,RCVD_IN_MSPIKE_H2,SPF_HELO_NONE,SPF_NONE,TXREP autolearn=ham autolearn_force=no version=3.4.6 X-Spam-Checker-Version: SpamAssassin 3.4.6 (2021-04-09) on server2.sourceware.org List-Id: On Thu, Nov 10, 2022 at 08:20:06PM +0100, Jakub Jelinek via Gcc-patches wrote: > 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). Ok, here is the patch rewritten in the foperator_div style, with special cases handled first and then the ordinary cases without problematic cases. I guess if/once we have a plugin testing infrastructure, we could compare the two versions of the patch, I think this one is more precise. And, admittedly there are many similar spots with the foperator_div case (but also with significant differences), so perhaps if foperator_{mult,div} inherit from some derived class from range_operator_float and that class would define various smaller helper static? methods, like this discussed in the PR - contains_zero_p, singleton_nan_p, zero_p, that + bool must_have_signbit_zero = false; + bool must_have_signbit_nonzero = 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; + } returned as -1/0/1 int, and those set result (based on the above value) to [+INF, +INF], [-INF, -INF] or [-INF, +INF] or [+0, +0], [-0, -0] or [-0, +0] or [+0, +INF], [-INF, -0] or [-INF, +INF] and the + 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]); + } block, it could be smaller and more readable. Thoughts? This has been just compile tested so far. 2022-11-11 Jakub Jelinek 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. (foperator_mult): New class. (floating_op_table::floating_op_table): Use foperator_mult for MULT_EXPR. --- gcc/range-op.h.jj 2022-11-11 08:15:20.952520590 +0100 +++ gcc/range-op.h 2022-11-11 08:48:27.649349048 +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-11 08:15:20.933520849 +0100 +++ gcc/range-op-float.cc 2022-11-11 09:39:14.950523368 +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,182 @@ class foperator_minus : public range_ope } } fop_minus; + +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 + { + 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)); + + maybe_nan = false; + // x * x never produces a new NAN and we only multiply the same + // values, so the 0 * INF problematic cases never appear there. + if (!is_square) + { + // [+-0, +-0] * [+INF,+INF] (or [-INF,-INF] or swapped is a known NAN. + 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; + return; + } + + // Otherwise, if one range includes zero and the other ends with +-INF, + // it is a maybe NAN. + if ((real_compare (LE_EXPR, &lh_lb, &dconst0) + && real_compare (GE_EXPR, &lh_ub, &dconst0) + && (real_isinf (&rh_lb) || real_isinf (&rh_ub))) + || (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; + + bool must_have_signbit_zero = false; + bool must_have_signbit_nonzero = 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 one of the ranges that includes INF is singleton + // and the other range includes zero, the resulting + // range is INF and NAN, because the 0 * INF boundary + // case will be NAN, but already nextafter (0, 1) * INF + // is INF. + if ((real_isinf (&lh_lb) + && real_isinf (&lh_ub, real_isneg (&lh_lb))) + || (real_isinf (&rh_lb) + && real_isinf (&rh_ub, real_isneg (&rh_lb)))) + { + // If all the boundary signs are the same, [+INF, +INF]. + if (must_have_signbit_zero) + ub = lb = dconstinf; + // If the two multiplicands have always different sign, + // [-INF, -INF]. + else if (must_have_signbit_nonzero) + ub = lb = dconstninf; + // Otherwise -> [-INF, +INF] (-INF or +INF). + else + { + lb = dconstninf; + ub = dconstinf; + } + return; + } + + // If one of the multiplicands must be zero, the resulting + // range is +-0 and NAN. + if ((real_iszero (&lh_lb) && real_iszero (&lh_ub)) + || (real_iszero (&rh_lb) && real_iszero (&rh_ub))) + { + 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; + } + + // Otherwise one of the multiplicands could be + // [0.0, nextafter (0.0, 1.0)] and the [DBL_MAX, INF] + // or similarly with different signs. 0.0 * DBL_MAX + // is still 0.0, nextafter (0.0, 1.0) * INF is still INF, + // so if the signs are always the same or always different, + // result is [+0.0, +INF] or [-INF, -0.0], otherwise VARYING. + 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-product. + frange_arithmetic (MULT_EXPR, type, cp[0], lh_lb, rh_lb, dconstninf); + frange_arithmetic (MULT_EXPR, type, cp[4], lh_lb, rh_lb, dconstinf); + 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_arithmetic (MULT_EXPR, type, cp[1], lh_lb, rh_ub, dconstninf); + frange_arithmetic (MULT_EXPR, type, cp[5], lh_lb, rh_ub, dconstinf); + frange_arithmetic (MULT_EXPR, type, cp[2], lh_ub, rh_lb, dconstninf); + frange_arithmetic (MULT_EXPR, type, cp[6], lh_ub, rh_lb, dconstinf); + } + frange_arithmetic (MULT_EXPR, type, cp[3], lh_ub, rh_ub, dconstninf); + frange_arithmetic (MULT_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]; + + } +} fop_mult; + // Instantiate a range_op_table for floating point operations. static floating_op_table global_floating_table; @@ -1942,6 +2120,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 Jakub