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

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

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

--- Comment #28 from Jakub Jelinek <jakub at gcc dot gnu.org> ---
On IRC we've discussed this and I believe a possible fix could be before we do:
    return float_binary_op_range_finish (minus.fold_range (r, type, lhs, op2),
                                         r, type, lhs);
etc. artificially extend the lhs range by 1ulp or .5ulp or whatever works in
each direction (on a copy) and then just let it do its job.

But I want to fixup above patch first, so that we have something to compare to.
From
--- gcc/range-op-float.cc.jj    2023-03-07 21:20:49.885225381 +0100
+++ gcc/range-op-float.cc       2023-03-08 09:13:10.063608296 +0100
@@ -2319,9 +2319,15 @@ float_range_extend (tree type, frange &r
        u = m;
     }
   gcc_checking_assert (l + 1 == u);
+fprintf (stderr, "--\n");
+SET_REAL_EXP (&t, l);
+real_arithmetic (&w, PLUS_EXPR, &v, &t);
+real_convert (&w, mode, &w);
+test (w, type, lhs, op2);
   SET_REAL_EXP (&t, u);
   real_arithmetic (&w, PLUS_EXPR, &v, &t);
   real_convert (&w, mode, &w);
+test (w, type, lhs, op2);
   REAL_VALUE_TYPE lastw = w;
   for (int i = SIGNIFICAND_BITS - 2; i >= SIGNIFICAND_BITS - 2 - p; i--)
     {
@@ -2338,7 +2344,10 @@ float_range_extend (tree type, frange &r
       lastw = w;
     }
   w = lastw;
+fprintf (stderr, "---\n");
+test (w, type, lhs, op2);
   frange_nextafter (mode, w, upper ? dconstninf : dconstinf);
+test (w, type, lhs, op2);
   goto update;
 }

@@ -2367,7 +2376,15 @@ public:
                            REAL_VALUE_TYPE res;
                            frange_arithmetic (PLUS_EXPR, type, res, r,
                                               op2.lower_bound (), dconstinf);
-                           return !real_less (&res, &lhs.lower_bound ());
+bool ret = !real_less (&res, &lhs.lower_bound ());
+char br[60], bop2[60], bres[60], blhs[60];
+real_to_hexadecimal (br, &r, sizeof (br), 0, 1);
+real_to_hexadecimal (bop2, &op2.lower_bound (), sizeof (bop2), 0, 1);
+real_to_hexadecimal (bres, &res, sizeof (bres), 0, 1);
+real_to_hexadecimal (blhs, &lhs.lower_bound (), sizeof (blhs), 0, 1);
+fprintf (stderr, "float_range_extend1 %s + %s = %s %s %s\n", br, bop2, bres,
ret ? ">=" : "<", blhs);
+return ret;
+//                         return !real_less (&res, &lhs.lower_bound ());
                          });
     float_range_extend (type, r, lhs, op2, true,
                        [] (REAL_VALUE_TYPE &r, tree type,
@@ -2376,7 +2393,15 @@ public:
                            REAL_VALUE_TYPE res;
                            frange_arithmetic (PLUS_EXPR, type, res, r,
                                               op2.upper_bound (), dconstninf);
-                           return !real_less (&lhs.upper_bound (), &res);
+bool ret = !real_less (&lhs.upper_bound (), &res);
+char br[60], bop2[60], bres[60], blhs[60];
+real_to_hexadecimal (br, &r, sizeof (br), 0, 1);
+real_to_hexadecimal (bop2, &op2.lower_bound (), sizeof (bop2), 0, 1);
+real_to_hexadecimal (bres, &res, sizeof (bres), 0, 1);
+real_to_hexadecimal (blhs, &lhs.lower_bound (), sizeof (blhs), 0, 1);
+fprintf (stderr, "float_range_extend2 %s + %s = %s %s %s\n", br, bop2, bres,
ret ? "<=" : ">", blhs);
+return ret;
+//                         return !real_less (&lhs.upper_bound (), &res);
                          });
     return true;
   }
debugging hacks seems it is the loop that tries to narrow the mantissa bits
that is wrong, as immediately before it the l and u values seem correct:
float_range_extend1 -0x0.8p-53 + 0x0.8p+1 = 0x0.8p+1 >= 0x0.8p+1
float_range_extend1 -0x0.8p-52 + 0x0.8p+1 = 0x0.fffffffffffff8p+0 < 0x0.8p+1
where the first is for l where test passes and u doesn't (and the right answer
here is -0x0.8p-53).