[Bug tree-optimization/109008] [13 Regression] Wrong code in scipy package since r13-3926-gd4c2f1d376da6f

["rguenth at gcc dot gnu.org" <gcc-bugzilla@gcc.gnu.org>]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=109008

--- Comment #13 from Richard Biener <rguenth at gcc dot gnu.org> ---
The question is what we want to do for GCC 13 - I suppose iterating would work
but it'll be slow (what's the range to binary search here?).  Doing the math
"right" probably differs for each reverse operation (but how many do we have?
I see plus, minus, mult and div only).

I'd say approximate but conservatively correct math would be prefered here.
Most conservative is probably to simply extend the input ranges by 1 ulp,
that should work for all operators but inevitably leads to less precise
answers.  As said, we could also resort to mpfr and widen the ranges by
0.5 ulp only, the question is whether that's less costly than iterating.
We could always at least verify if the forward operation with the result
yields something that contains the old LHS range though (maybe that's a
good idea anyway when checking is enabled).

But the real challenge of course is to somehow have test coverage for all this.

With something like

diff --git a/gcc/range-op-float.cc b/gcc/range-op-float.cc
index ff42b95de4f..e142b83c047 100644
--- a/gcc/range-op-float.cc
+++ b/gcc/range-op-float.cc
@@ -2214,7 +2214,29 @@ 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),
+    /* ???  some first-class "widening" CTOR might be nicer, maybe
+       add some static function?  */
+    frange wlhs (lhs);
+    if (!lhs.known_isnan ())
+      {
+       REAL_VALUE_TYPE lhsl = lhs.lower_bound ();
+       frange_nextafter (TYPE_MODE (type), lhsl, dconstninf);
+       REAL_VALUE_TYPE lhsu = lhs.upper_bound ();
+       frange_nextafter (TYPE_MODE (type), lhsu, dconstinf);
+       wlhs.set (type, lhsl, lhsu);
+       // no way to copy NaN state?
+      }
+    frange wop2 (op2);
+    if (!op2.known_isnan ())
+      {
+       REAL_VALUE_TYPE op2l = op2.lower_bound ();
+       frange_nextafter (TYPE_MODE (type), op2l, dconstninf);
+       REAL_VALUE_TYPE op2u = op2.upper_bound ();
+       frange_nextafter (TYPE_MODE (type), op2u, dconstinf);
+       wop2.set (type, op2l, op2u);
+       // no way to copy NaN state?
+      }
+    return float_binary_op_range_finish (minus.fold_range (r, type, wlhs,
wop2),
                                         r, type, lhs);
   }
   virtual bool op2_range (frange &r, tree type,

the comment#2 testcase shows in EVRP

Imports: eps_2(D)
Exports: eps_2(D)  d_3
         d_3 : eps_2(D)(I)
eps_2(D)        [frange] double VARYING +-NAN
    <bb 2> :
    d_3 = eps_2(D) + 1.0e+0;
    if (d_3 == 1.0e+0)
      goto <bb 3>; [INV]
    else
      goto <bb 4>; [INV]

2->3  (T) eps_2(D) :    [frange] double
[-3.3306690738754696212708950042724609375e-16 (-0x0.cp-51),
3.3306690738754696212708950042724609375e-16 (0x0.cp-51)]
2->3  (T) d_3 :         [frange] double [1.0e+0 (0x0.8p+1), 1.0e+0 (0x0.8p+1)]

that's at least no longer incorrect.

As written in the comment above mangling a range like this is currently
a bit awkward, but maybe a static method in frange like
frange frange::1ulp_wider (const frange &) would be possible.

If we ever get range ops for things like exp() things will get more
interesting, the input range widening by n * ulp only works for linear
ops I think.