From patchwork Thu Mar 9 08:29:02 2023 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Jakub Jelinek X-Patchwork-Id: 66155 Return-Path: X-Original-To: patchwork@sourceware.org Delivered-To: patchwork@sourceware.org Received: from server2.sourceware.org (localhost [IPv6:::1]) by sourceware.org (Postfix) with ESMTP id DE5CA385802F for ; Thu, 9 Mar 2023 08:29:42 +0000 (GMT) DKIM-Filter: OpenDKIM Filter v2.11.0 sourceware.org DE5CA385802F DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gcc.gnu.org; s=default; t=1678350582; bh=1eyGadaxc9/7KTO/C9iLsTqUl5U2Kao4Qbu+6ali6cA=; h=Date:To:Cc:Subject:List-Id:List-Unsubscribe:List-Archive: List-Post:List-Help:List-Subscribe:From:Reply-To:From; b=w3muRQIUOzOoReUDHTuVYfeTg28Tufz9VOGCAV2hw1McOxFFDufDCTDu3jX9BhTEh cSO58VW/KMRzAwc33l/lkipNjv6fL7br56/vE1lhEeiq8k/N0wKFtuGkOLlfxW0+3T RmCW2rKlP585c0ab4aCv9XEB5tLEL3r9nBPXlbmU= X-Original-To: gcc-patches@gcc.gnu.org Delivered-To: gcc-patches@gcc.gnu.org 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 78040385B51B for ; Thu, 9 Mar 2023 08:29:08 +0000 (GMT) DMARC-Filter: OpenDMARC Filter v1.4.2 sourceware.org 78040385B51B 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-620-J0D03j5EPu2g_j_0WVWh5g-1; Thu, 09 Mar 2023 03:29:07 -0500 X-MC-Unique: J0D03j5EPu2g_j_0WVWh5g-1 Received: from smtp.corp.redhat.com (int-mx04.intmail.prod.int.rdu2.redhat.com [10.11.54.4]) (using TLSv1.2 with cipher AECDH-AES256-SHA (256/256 bits)) (No client certificate requested) by mimecast-mx02.redhat.com (Postfix) with ESMTPS id ACF86185A78B; Thu, 9 Mar 2023 08:29:06 +0000 (UTC) Received: from tucnak.zalov.cz (unknown [10.39.192.16]) by smtp.corp.redhat.com (Postfix) with ESMTPS id 6E0BE2026D4B; Thu, 9 Mar 2023 08:29:06 +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 3298T3PP695675 (version=TLSv1.3 cipher=TLS_AES_256_GCM_SHA384 bits=256 verify=NOT); Thu, 9 Mar 2023 09:29:04 +0100 Received: (from jakub@localhost) by tucnak.zalov.cz (8.17.1/8.17.1/Submit) id 3298T38w695674; Thu, 9 Mar 2023 09:29:03 +0100 Date: Thu, 9 Mar 2023 09:29:02 +0100 To: Richard Biener , Aldy Hernandez Cc: gcc-patches@gcc.gnu.org Subject: [PATCH] range-op-float: Fix up reverse binary operations [PR109008] Message-ID: MIME-Version: 1.0 X-Scanned-By: MIMEDefang 3.1 on 10.11.54.4 X-Mimecast-Spam-Score: 0 X-Mimecast-Originator: redhat.com Content-Disposition: inline X-Spam-Status: No, score=-3.4 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 X-BeenThere: gcc-patches@gcc.gnu.org X-Mailman-Version: 2.1.29 Precedence: list List-Id: Gcc-patches mailing list List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-Patchwork-Original-From: Jakub Jelinek via Gcc-patches From: Jakub Jelinek Reply-To: Jakub Jelinek Errors-To: gcc-patches-bounces+patchwork=sourceware.org@gcc.gnu.org Sender: "Gcc-patches" Hi! The following testcase is reduced from miscompilation of scipy package. If we have say lhs = [1., 1.] - [1., 1.] and want to compute the range of lhs from it, we correctly determine it is [0., 0.] (if computations are exact, we generally don't try to round them further in frange_arithmetic). In the testcase it is about a reverse operation, [1., 1.] = op1 + [1., 1.] and we want to compute range of op1 from that. Right now we just perform the inverse operation (there are some corner cases about NaN and infinities handling) and so arrive to range [0., 0.] as well, and because it is a singleton, optimize return eps; to return 0. That is incorrect though, for the reverse ops we need to take into account also rounding, the right exact range is [-0x1.0p-54, 0x1.0p-53] in this case when rounding to nearest, i.e. all numbers which added to 1. with round to nearest still produce 1. The problem isn't solely on singleton ranges, and isn't solely on results around zero. We basically need to consider also values where the result is up to 0.5ulp away from the lhs range boundaries in each direction. The following patch fixes it by extending the lhs range for the reverse operations by 1ulp in each direction. The PR contains a pseudo-random test generator I've used to generate 300000 tests of + and - and then used the same test with * and / instead of + and - together with a hack to print the discovered ranges by the patch in a form that another test could then verify the range is conservatively correct and how far it is from a minimal range. I believe the results are good enough for now, though plan to look incrementally into trying to do something better on the -XXX_MAX or XXX_MAX boundaries (where I think frange_nextafter will use -inf or +inf) and also try to increase the range just by 0.5ulp rather than 1ulp if !flag_rounding_math. But dunno if either of those will be doable and will pass the testing, so I think it is worth committing this fix first. Bootstrapped/regtested on x86_64-linux and i686-linux, ok for trunk? 2023-03-09 Jakub Jelinek Richard Biener PR tree-optimization/109008 * range-op-float.cc (float_widen_lhs_range): New function. (foperator_plus::op1_range, foperator_minus::op1_range, foperator_minus::op2_range, foperator_mult::op1_range, foperator_div::op1_range, foperator_div::op2_range): Use it. * gcc.c-torture/execute/ieee/pr109008.c: New test. Jakub --- gcc/range-op-float.cc.jj 2023-03-08 12:33:44.641043477 +0100 +++ gcc/range-op-float.cc 2023-03-08 13:13:09.015341002 +0100 @@ -2199,6 +2199,33 @@ zero_to_inf_range (REAL_VALUE_TYPE &lb, } } +/* Extend the LHS range by 1ulp in each direction. For op1_range + or op2_range of binary operations just computing the inverse + operation on ranges isn't sufficient. Consider e.g. + [1., 1.] = op1 + [1., 1.]. op1's range is not [0., 0.], but + [-0x1.0p-54, 0x1.0p-53] (when not -frounding-math), any value for + which adding 1. to it results in 1. after rounding to nearest. + So, for op1_range/op2_range extend the lhs range by 1ulp in each + direction. See PR109008 for more details. */ + +static frange +float_widen_lhs_range (tree type, const frange &lhs) +{ + frange ret = lhs; + if (lhs.known_isnan ()) + return ret; + REAL_VALUE_TYPE lb = lhs.lower_bound (); + REAL_VALUE_TYPE ub = lhs.upper_bound (); + if (real_isfinite (&lb)) + frange_nextafter (TYPE_MODE (type), lb, dconstninf); + if (real_isfinite (&ub)) + frange_nextafter (TYPE_MODE (type), ub, dconstinf); + ret.set (type, lb, ub); + ret.clear_nan (); + ret.union_ (lhs); + return ret; +} + class foperator_plus : public range_operator_float { using range_operator_float::op1_range; @@ -2214,8 +2241,9 @@ public: range_op_handler minus (MINUS_EXPR, type); if (!minus) return false; - return float_binary_op_range_finish (minus.fold_range (r, type, lhs, op2), - r, type, lhs); + frange wlhs = float_widen_lhs_range (type, lhs); + return float_binary_op_range_finish (minus.fold_range (r, type, wlhs, op2), + r, type, wlhs); } virtual bool op2_range (frange &r, tree type, const frange &lhs, @@ -2260,9 +2288,10 @@ public: { if (lhs.undefined_p ()) return false; - return float_binary_op_range_finish (fop_plus.fold_range (r, type, lhs, + frange wlhs = float_widen_lhs_range (type, lhs); + return float_binary_op_range_finish (fop_plus.fold_range (r, type, wlhs, op2), - r, type, lhs); + r, type, wlhs); } virtual bool op2_range (frange &r, tree type, const frange &lhs, @@ -2271,8 +2300,9 @@ public: { if (lhs.undefined_p ()) return false; - return float_binary_op_range_finish (fold_range (r, type, op1, lhs), - r, type, lhs); + frange wlhs = float_widen_lhs_range (type, lhs); + return float_binary_op_range_finish (fold_range (r, type, op1, wlhs), + r, type, wlhs); } private: void rv_fold (REAL_VALUE_TYPE &lb, REAL_VALUE_TYPE &ub, bool &maybe_nan, @@ -2338,13 +2368,14 @@ public: range_op_handler rdiv (RDIV_EXPR, type); if (!rdiv) return false; - bool ret = rdiv.fold_range (r, type, lhs, op2); + frange wlhs = float_widen_lhs_range (type, lhs); + bool ret = rdiv.fold_range (r, type, wlhs, op2); if (ret == false) return false; - if (lhs.known_isnan () || op2.known_isnan () || op2.undefined_p ()) - return float_binary_op_range_finish (ret, r, type, lhs); - const REAL_VALUE_TYPE &lhs_lb = lhs.lower_bound (); - const REAL_VALUE_TYPE &lhs_ub = lhs.upper_bound (); + if (wlhs.known_isnan () || op2.known_isnan () || op2.undefined_p ()) + return float_binary_op_range_finish (ret, r, type, wlhs); + const REAL_VALUE_TYPE &lhs_lb = wlhs.lower_bound (); + const REAL_VALUE_TYPE &lhs_ub = wlhs.upper_bound (); const REAL_VALUE_TYPE &op2_lb = op2.lower_bound (); const REAL_VALUE_TYPE &op2_ub = op2.upper_bound (); if ((contains_zero_p (lhs_lb, lhs_ub) && contains_zero_p (op2_lb, op2_ub)) @@ -2363,7 +2394,7 @@ public: // or if lhs must be zero and op2 doesn't include zero, it would be // UNDEFINED, while rdiv.fold_range computes a zero or singleton INF // range. Those are supersets of UNDEFINED, so let's keep that way. - return float_binary_op_range_finish (ret, r, type, lhs); + return float_binary_op_range_finish (ret, r, type, wlhs); } virtual bool op2_range (frange &r, tree type, const frange &lhs, @@ -2490,13 +2521,14 @@ public: { if (lhs.undefined_p ()) return false; - bool ret = fop_mult.fold_range (r, type, lhs, op2); + frange wlhs = float_widen_lhs_range (type, lhs); + bool ret = fop_mult.fold_range (r, type, wlhs, op2); if (!ret) return ret; - if (lhs.known_isnan () || op2.known_isnan () || op2.undefined_p ()) - return float_binary_op_range_finish (ret, r, type, lhs); - const REAL_VALUE_TYPE &lhs_lb = lhs.lower_bound (); - const REAL_VALUE_TYPE &lhs_ub = lhs.upper_bound (); + if (wlhs.known_isnan () || op2.known_isnan () || op2.undefined_p ()) + return float_binary_op_range_finish (ret, r, type, wlhs); + const REAL_VALUE_TYPE &lhs_lb = wlhs.lower_bound (); + const REAL_VALUE_TYPE &lhs_ub = wlhs.upper_bound (); const REAL_VALUE_TYPE &op2_lb = op2.lower_bound (); const REAL_VALUE_TYPE &op2_ub = op2.upper_bound (); if ((contains_zero_p (lhs_lb, lhs_ub) @@ -2512,7 +2544,7 @@ public: zero_to_inf_range (lb, ub, signbit_known); r.set (type, lb, ub); } - return float_binary_op_range_finish (ret, r, type, lhs); + return float_binary_op_range_finish (ret, r, type, wlhs); } virtual bool op2_range (frange &r, tree type, const frange &lhs, @@ -2521,13 +2553,14 @@ public: { if (lhs.undefined_p ()) return false; - bool ret = fold_range (r, type, op1, lhs); + frange wlhs = float_widen_lhs_range (type, lhs); + bool ret = fold_range (r, type, op1, wlhs); if (!ret) return ret; - if (lhs.known_isnan () || op1.known_isnan () || op1.undefined_p ()) - return float_binary_op_range_finish (ret, r, type, lhs, true); - const REAL_VALUE_TYPE &lhs_lb = lhs.lower_bound (); - const REAL_VALUE_TYPE &lhs_ub = lhs.upper_bound (); + if (wlhs.known_isnan () || op1.known_isnan () || op1.undefined_p ()) + return float_binary_op_range_finish (ret, r, type, wlhs, true); + const REAL_VALUE_TYPE &lhs_lb = wlhs.lower_bound (); + const REAL_VALUE_TYPE &lhs_ub = wlhs.upper_bound (); const REAL_VALUE_TYPE &op1_lb = op1.lower_bound (); const REAL_VALUE_TYPE &op1_ub = op1.upper_bound (); if ((contains_zero_p (lhs_lb, lhs_ub) && contains_zero_p (op1_lb, op1_ub)) @@ -2542,7 +2575,7 @@ public: zero_to_inf_range (lb, ub, signbit_known); r.set (type, lb, ub); } - return float_binary_op_range_finish (ret, r, type, lhs, true); + return float_binary_op_range_finish (ret, r, type, wlhs, true); } private: void rv_fold (REAL_VALUE_TYPE &lb, REAL_VALUE_TYPE &ub, bool &maybe_nan, --- gcc/testsuite/gcc.c-torture/execute/ieee/pr109008.c.jj 2023-03-08 21:30:19.158618157 +0100 +++ gcc/testsuite/gcc.c-torture/execute/ieee/pr109008.c 2023-03-08 21:29:49.899039474 +0100 @@ -0,0 +1,18 @@ +/* PR tree-optimization/109008 */ + +__attribute__((noipa)) double +foo (double eps) +{ + double d = 1. + eps; + if (d == 1.) + return eps; + return 0.0; +} + +int +main () +{ + if (foo (__DBL_EPSILON__ / 8.0) == 0.0) + __builtin_abort (); + return 0; +}