Re: [PATCH] [PR24021] Implement PLUS_EXPR range-op entry for floats.

[Aldy Hernandez <aldyh@redhat.com>]

PING

On Mon, Oct 17, 2022 at 8:21 AM Aldy Hernandez <aldyh@redhat.com> wrote:
>
> On Thu, Oct 13, 2022 at 7:57 PM Jakub Jelinek <jakub@redhat.com> wrote:
> >
> > On Thu, Oct 13, 2022 at 02:36:49PM +0200, Aldy Hernandez wrote:
> > > +// Like real_arithmetic, but round the result to INF if the operation
> > > +// produced inexact results.
> > > +//
> > > +// ?? There is still one problematic case, i387.  With
> > > +// -fexcess-precision=standard we perform most SF/DFmode arithmetic in
> > > +// XFmode (long_double_type_node), so that case is OK.  But without
> > > +// -mfpmath=sse, all the SF/DFmode computations are in XFmode
> > > +// precision (64-bit mantissa) and only occassionally rounded to
> > > +// SF/DFmode (when storing into memory from the 387 stack).  Maybe
> > > +// this is ok as well though it is just occassionally more precise. ??
> > > +
> > > +static void
> > > +frange_arithmetic (enum tree_code code, tree type,
> > > +                REAL_VALUE_TYPE &result,
> > > +                const REAL_VALUE_TYPE &op1,
> > > +                const REAL_VALUE_TYPE &op2,
> > > +                const REAL_VALUE_TYPE &inf)
> > > +{
> > > +  REAL_VALUE_TYPE value;
> > > +  enum machine_mode mode = TYPE_MODE (type);
> > > +  bool mode_composite = MODE_COMPOSITE_P (mode);
> > > +
> > > +  bool inexact = real_arithmetic (&value, code, &op1, &op2);
> > > +  real_convert (&result, mode, &value);
> > > +
> > > +  // If real_convert above has rounded an inexact value to towards
> > > +  // inf, we can keep the result as is, otherwise we'll adjust by 1 ulp
> > > +  // later (real_nextafter).
> > > +  bool rounding = (flag_rounding_math
> > > +                && (real_isneg (&inf)
> > > +                    ? real_less (&result, &value)
> > > +                    : !real_less (&value, &result)));
> >
> > I thought the agreement during Cauldron was that we'd do this always,
> > regardless of flag_rounding_math.
> > Because excess precision (the fast one like on ia32 or -mfpmath=387 on
> > x86_64), or -frounding-math, or FMA contraction can all increase precision
> > and worst case it all behaves like -frounding-math for the ranges.
> >
> > So, perhaps use:
> >   if ((mode_composite || (real_isneg (&inf) ? real_less (&result, &value)
> >                                             : !real_less (&value, &result))
> >       && (inexact || !real_identical (&result, &value))))
>
> Done.
>
> > ?
> > No need to do the real_isneg/real_less stuff for mode_composite, then
> > we do it always for inexacts, but otherwise we check if the rounding
> > performed by real.cc has been in the conservative direction (for upper
> > bound to +inf, for lower bound to -inf), if yes, we don't need to do
> > anything, if yes, we frange_nextafter.
> >
> > As discussed, for mode_composite, I think we want to do the extra
> > stuff for inexact denormals and otherwise do the nextafter unconditionally,
> > because our internal mode_composite representation isn't precise enough.
> >
> > > +  // Be extra careful if there may be discrepancies between the
> > > +  // compile and runtime results.
> > > +  if ((rounding || mode_composite)
> > > +      && (inexact || !real_identical (&result, &value)))
> > > +    {
> > > +      if (mode_composite)
> > > +     {
> > > +       bool denormal = (result.sig[SIGSZ-1] & SIG_MSB) == 0;
> >
> > Use real_isdenormal here?
>
> Done.
>
> > Though, real_iszero needs the same thing.
>
> So... real_isdenormal() || real_iszero() as in the attached patch?
>
> >
> > > +       if (denormal)
> > > +         {
> > > +           REAL_VALUE_TYPE tmp;
> >
> > And explain here why is this, that IBM extended denormals have just
> > DFmode precision.
>
> Done.
>
> > Though, now that I think about it, while this is correct for denormals,
> >
> > > +           real_convert (&tmp, DFmode, &value);
> > > +           frange_nextafter (DFmode, tmp, inf);
> > > +           real_convert (&result, mode, &tmp);
> > > +         }
> >
> > there are also the cases where the higher double exponent is in the
> > [__DBL_MIN_EXP__, __LDBL_MIN_EXP__] aka [-1021, -968] or so.
> > https://en.wikipedia.org/wiki/Double-precision_floating-point_format
> > If the upper double is denormal in the DFmode sense, so smaller absolute
> > value than __DBL_MIN__, then doing nextafter in DFmode is the right thing to
> > do, the lower double must be always +/- zero.
> > Now, if the result is __DBL_MIN__, the upper double is already normalized
> > but we can add __DBL_DENORM_MIN__ to it, which will make the number have
> > 54-bit precision.
> > If the result is __DBL_MIN__ * 2, we can again add __DBL_DENORM_MIN__
> > and make it 55-bit precision.  Etc. until we reach __DBL_MIN__ * 2e53
> > where it acts like fully normalized 106-bit precision number.
> > I must say I'm not really sure what real_nextafter is doing in those cases,
> > I'm afraid it doesn't handle it correctly but the only other use
> > of real_nextafter is guarded with:
> >   /* Don't handle composite modes, nor decimal, nor modes without
> >      inf or denorm at least for now.  */
> >   if (format->pnan < format->p
> >       || format->b == 10
> >       || !format->has_inf
> >       || !format->has_denorm)
> >     return false;
>
> Dunno.  Is there a conservative thing we can do for mode_composites
> that aren't denormal or zero?
>
> How does this look?
> Aldy