https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107569
--- Comment #18 from Jakub Jelinek <jakub at gcc dot gnu.org> --- Ok, just so that I don't just kibbitz/review frange stuff but also try to do something, here is my so far untested attempt at normal multiplication fold_range (not the x * x stuff discussed elsewhere): --- range-op-float.cc.jj1 2022-11-09 18:00:28.612247613 +0100 +++ range-op-float.cc 2022-11-09 19:06:11.075716000 +0100 @@ -1908,6 +1908,88 @@ 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. */ + +static void +frange_mult (tree type, REAL_VALUE_TYPE &result_lb, REAL_VALUE_TYPE &result_ub, + const REAL_VALUE_TYPE &op1, const REAL_VALUE_TYPE &op2) +{ + if (real_iszero (&op1) && real_isinf (&op2)) + { + result_lb = op1; + if (real_isneg (&op2)) + real_value_negate (&result_lb); + result_ub = result_lb; + } + else if (real_isinf (&op1) && real_iszero (&op2)) + { + result_lb = op2; + if (real_isneg (&op1)) + real_value_negate (&result_lb); + result_ub = result_lb; + } + else + { + frange_arithmetic (MULT_EXPR, type, result_lb, op1, op2, dconstninf); + frange_arithmetic (MULT_EXPR, type, result_ub, op1, op2, dconstinf); + } +} + +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) const final override + { + REAL_VALUE_TYPE cp[8]; + // Do a cross-product. + frange_mult (type, cp[0], cp[4], lh_lb, rh_lb); + frange_mult (type, cp[1], cp[5], lh_lb, rh_ub); + frange_mult (type, cp[2], cp[6], lh_ub, rh_lb); + frange_mult (type, cp[3], cp[7], lh_ub, rh_ub); + for (int i = 1; i < 3; ++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]; + + // [+-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, NULL, 0, TYPE_MODE (type)); + ub = lb; + maybe_nan = true; + } + // Otherwise, if one range includes zero and the other ends with +-INF, + // it is a maybe NaN. + else if (real_compare (LE_EXPR, &lh_lb, &dconst0) + && real_compare (GE_EXPR, &lh_ub, &dconst0) + && (real_isinf (&rh_lb) || real_isinf (&rh_ub))) + maybe_nan = true; + else if (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; + else + maybe_nan = false; + } +} fop_mult; + // Instantiate a range_op_table for floating point operations. static floating_op_table global_floating_table; @@ -1942,6 +2024,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 For float foo (float x, float y) { if (x < -17.0f || x > 19.5f || y < -25.0f || y > 15.5f) return 0.0; return x * y; } float bar (float x, float y) { if (x > -17.0f && x < 19.5f && y > -25.0f && y < 15.5f) return x * y; return 0.0; } the product range is [frange] float [-4.875e+2 (-0x0.f3cp+9), 4.25e+2 (0x0.d48p+9)] +-NAN in the first case and [frange] float [-4.875e+2 (-0x0.f3cp+9), 4.25e+2 (0x0.d48p+9)] in the latter (where NAN can't appear). And while +-0.0 is in the range of even both operands, neither range includes any infinities.