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)))) ? 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? Though, real_iszero needs the same thing. > + if (denormal) > + { > + REAL_VALUE_TYPE tmp; And explain here why is this, that IBM extended denormals have just DFmode precision. 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; so it isn't that big deal except for ranges. Jakub