https://gcc.gnu.org/bugzilla/show_bug.cgi?id=127119
Bug ID: 127119
Summary: add vectorization for 'round' with flag_trapping_math
Product: gcc
Version: 17.0
Status: UNCONFIRMED
Keywords: missed-optimization
Severity: enhancement
Priority: P3
Component: target
Assignee: unassigned at gcc dot gnu.org
Reporter: mkretz at gcc dot gnu.org
CC: crazylht at gmail dot com, hongyuw at gcc dot gnu.org
Target Milestone: ---
Target: x86_64-*-*
As discussed in PR127055, round can be vectorized even with -ftrapping-math.
However in PR127055#c19 I proposed an incorrect solution:
abs_x = fabs(x);
t_abs = trunc(abs_x);
r_abs = t_abs + (abs_x - t_abs >= .5 ? 1 : 0);
return or(xor(abs_x, x), r_abs);
This is wrong at x=inf, because abs_x - t_abs raises FE_INVALID. (I fixed my
unit tests to find this now.)
I then rewrote the 'abs_x - t_abs >= .5' test using integer operations. The
resulting function passes all tests (all fp exceptions match exactly what
scalar std::round does). Even more interesting, according to llvm-mca:
-march= / --mcpu= | __round | current -Ofast
--------------------|-----------|---------------------------
skylake / skylake | 19 cycles | 14 cycles
znver5 / znver5 | 11 cycles | 14 cycles
x86-64-v3 / znver5 | 7 cycles | 8 cycles
This is the C++ implementation (using libstdc++-v3/include/bits/vec_ops.h
helpers):
template <__vec_builtin _TV>
_TV
__round(_TV __x)
{
using _Tp = __vec_value_type<_TV>;
using _Up = _UInt<sizeof(_Tp)>;
using _Ip = __integer_from<sizeof(_Tp)>;
const auto __abs_x = __fabs<_Traits>(__x);
const auto __t_abs = __trunc<_Traits>(__abs_x);
const auto __int_abs_x = __vec_bit_cast<_Ip>(__abs_x);
const auto __int_t_abs = __vec_bit_cast<_Ip>(__t_abs);
// a) for |x| without fractional mantissa bits: __diff = 0
// b) for |x| >= 1: XOR cancels exponent bits and mantissa bits
signifying values >= 1
// If the MSB is at the .5 position we need to add 1 to
the result
// c) for |x| < 1: trunc(|x|) is 0 => so __diff = bit-pattern of |x|
// => __diff needs to be compared against
// a) any positive non-zero integer
// b) integer with 1 bit at .5 position
// c) bit-pattern of .5
const auto __diff = __int_abs_x ^ __int_t_abs;
constexpr int __mant_width = numeric_limits<_Tp>::digits - 1;
constexpr _Ip __one_bits = __builtin_bit_cast(_Ip, _Tp(1));
constexpr _Ip __one_half_bits = __builtin_bit_cast(_Ip, _Tp(.5));
constexpr _Ip __bias = __one_bits >> __mant_width;
const auto __biased_exp = __int_abs_x >> __mant_width;
const auto __exponent =
__vec_bit_cast<_Ip>(__vec_bit_cast<_Up>(__biased_exp - __bias) %
numeric_limits<_Up>::digits);
const auto __threshold
= ((__biased_exp > (__mant_width + __bias)) | (__int_abs_x <
__one_bits))
? __one_half_bits // a) and c)
: (_Ip(1) << (__mant_width - 1)) >> __exponent; // b)
const _TV __sign_bit = __vec_xor(__abs_x, __x);
const _TV __r_abs = __t_abs + __builtin_bit_cast(
_TV, __diff < __threshold ? 0 :
__one_bits);
return __vec_or(__sign_bit, __r_abs);
}
I replaced fp subtract with XOR, doing integer compare against a value computed
in the "free" 8 cycles __trunc needs to compute __t_abs.
The main issue is that `>> __exponent` is a vector shift, which is only really
efficient with AVX2. However, improving vector shifts without AVX2 is again a
separate issue to be resolved.
Now, all that said, the resulting number of instructions is likely too much for
inlining. So I assume we rather want to have such a solution in the library,
right?
