From 716c2d08893110f186f0fb94849524b1acf9dd02 Mon Sep 17 00:00:00 2001 From: Jakub Jelinek Date: Thu, 8 Dec 2022 10:41:49 +0100 Subject: [PATCH] range-op-float: frange_arithmetic tweaks for MODE_COMPOSITE_P As mentioned in PR107967, ibm-ldouble-format documents that +- has 1ulp accuracy, * 2ulps and / 3ulps. So, even if the result is exact, we need to widen the range a little bit. The following patch does that. I just wonder what it means for reverse division (the op1_range case), which we implement through multiplication, when division has 3ulps error and multiplication just 2ulps. In any case, this format is a mess and for non-default rounding modes can't be trusted at all, instead of +inf or something close to it it happily computes -inf. 2022-12-08 Jakub Jelinek * range-op-float.cc (frange_nextafter): For MODE_COMPOSITE_P from denormal or zero, use real_nextafter on DFmode with conversions around it. (frange_arithmetic): For mode_composite, on top of rounding in the right direction accept extra 1ulp error for PLUS/MINUS_EXPR, extra 2ulps error for MULT_EXPR and extra 3ulps error for RDIV_EXPR. --- gcc/range-op-float.cc | 62 ++++++++++++++++++++++++++++++++++++++------------- 1 file changed, 46 insertions(+), 16 deletions(-) diff --git a/gcc/range-op-float.cc b/gcc/range-op-float.cc index 2c6026d..929ff9c 100644 --- a/gcc/range-op-float.cc +++ b/gcc/range-op-float.cc @@ -254,10 +254,21 @@ frange_nextafter (enum machine_mode mode, REAL_VALUE_TYPE &value, const REAL_VALUE_TYPE &inf) { - const real_format *fmt = REAL_MODE_FORMAT (mode); - REAL_VALUE_TYPE tmp; - real_nextafter (&tmp, fmt, &value, &inf); - value = tmp; + if (MODE_COMPOSITE_P (mode) + && (real_isdenormal (&value, mode) || real_iszero (&value))) + { + // IBM extended denormals only have DFmode precision. + REAL_VALUE_TYPE tmp, tmp2; + real_convert (&tmp2, DFmode, &value); + real_nextafter (&tmp, REAL_MODE_FORMAT (DFmode), &tmp2, &inf); + real_convert (&value, mode, &tmp); + } + else + { + REAL_VALUE_TYPE tmp; + real_nextafter (&tmp, REAL_MODE_FORMAT (mode), &value, &inf); + value = tmp; + } } // Like real_arithmetic, but round the result to INF if the operation @@ -324,21 +335,40 @@ frange_arithmetic (enum tree_code code, tree type, } if (round && (inexact || !real_identical (&result, &value))) { - if (mode_composite) + if (mode_composite + && (real_isdenormal (&result, mode) || real_iszero (&result))) { - if (real_isdenormal (&result, mode) - || real_iszero (&result)) - { - // IBM extended denormals only have DFmode precision. - REAL_VALUE_TYPE tmp; - real_convert (&tmp, DFmode, &value); - frange_nextafter (DFmode, tmp, inf); - real_convert (&result, mode, &tmp); - return; - } + // IBM extended denormals only have DFmode precision. + REAL_VALUE_TYPE tmp, tmp2; + real_convert (&tmp2, DFmode, &value); + real_nextafter (&tmp, REAL_MODE_FORMAT (DFmode), &tmp2, &inf); + real_convert (&result, mode, &tmp); } - frange_nextafter (mode, result, inf); + else + frange_nextafter (mode, result, inf); } + if (mode_composite) + switch (code) + { + case PLUS_EXPR: + case MINUS_EXPR: + // ibm-ldouble-format documents 1ulp for + and -. + frange_nextafter (mode, result, inf); + break; + case MULT_EXPR: + // ibm-ldouble-format documents 2ulps for *. + frange_nextafter (mode, result, inf); + frange_nextafter (mode, result, inf); + break; + case RDIV_EXPR: + // ibm-ldouble-format documents 3ulps for /. + frange_nextafter (mode, result, inf); + frange_nextafter (mode, result, inf); + frange_nextafter (mode, result, inf); + break; + default: + break; + } } // Crop R to [-INF, MAX] where MAX is the maximum representable number -- 2.7.4