From 871b0cec0f24b2a21bab1fc45ae8031cbadd7c0e Mon Sep 17 00:00:00 2001 From: Vladimir Saraikin Date: Mon, 13 Jul 2026 19:11:18 +0200 Subject: [PATCH] Fix saturating subtraction when the subtrahend is the type minimum ssub was implemented as sadd(x, -rhs), but -rhs is not representable when rhs == numeric_limits::min(), so e.g. ssub(1, INT_MIN) wrapped to a large negative value instead of saturating to INT_MAX. Compute the saturating difference directly: - scalar overload and the common batch kernel (used by SSE2 for 32/64-bit integers, which have no native saturating-subtract intrinsic); - AVX and AVX512F integer ssub, which had their own copies of the sadd(x, -rhs) pattern, now mirror their local sadd implementations. For 8/16-bit integers the saturating-sub intrinsic requires AVX512BW; in the AVX512F fallback, delegate those widths to AVX2 on each 256-bit half. --- .../arch/common/xsimd_common_arithmetic.hpp | 7 ++++++- include/xsimd/arch/xsimd_avx.hpp | 5 ++++- include/xsimd/arch/xsimd_avx512f.hpp | 15 +++++++++++++-- include/xsimd/arch/xsimd_scalar.hpp | 16 +++++++++++++++- test/test_xsimd_api.cpp | 19 +++++++++++++++++++ 5 files changed, 57 insertions(+), 5 deletions(-) diff --git a/include/xsimd/arch/common/xsimd_common_arithmetic.hpp b/include/xsimd/arch/common/xsimd_common_arithmetic.hpp index 541503ea3..b0ce69d2c 100644 --- a/include/xsimd/arch/common/xsimd_common_arithmetic.hpp +++ b/include/xsimd/arch/common/xsimd_common_arithmetic.hpp @@ -398,7 +398,12 @@ namespace xsimd { if (std::is_signed::value) { - return sadd(self, -other); + // Saturating self - other, mirroring the signed sadd above. + // sadd(self, -other) is wrong when other == numeric_limits::min(), + // since -other is not representable. + auto self_underflow_branch = max(std::numeric_limits::min() + other, self); + auto self_overflow_branch = min(std::numeric_limits::max() + other, self); + return select(other >= 0, self_underflow_branch, self_overflow_branch) - other; } else { diff --git a/include/xsimd/arch/xsimd_avx.hpp b/include/xsimd/arch/xsimd_avx.hpp index e048e4d37..ff95ca521 100644 --- a/include/xsimd/arch/xsimd_avx.hpp +++ b/include/xsimd/arch/xsimd_avx.hpp @@ -1711,7 +1711,10 @@ namespace xsimd { if (std::is_signed::value) { - return sadd(self, -other); + auto mask = (other >> (8 * sizeof(T) - 1)); + auto self_overflow_branch = min(std::numeric_limits::max() + other, self); + auto self_underflow_branch = max(std::numeric_limits::min() + other, self); + return select(batch_bool(mask.data), self_overflow_branch, self_underflow_branch) - other; } else { diff --git a/include/xsimd/arch/xsimd_avx512f.hpp b/include/xsimd/arch/xsimd_avx512f.hpp index b8a56c7ba..5500f39c7 100644 --- a/include/xsimd/arch/xsimd_avx512f.hpp +++ b/include/xsimd/arch/xsimd_avx512f.hpp @@ -2419,9 +2419,20 @@ namespace xsimd template ::value>> XSIMD_INLINE batch ssub(batch const& self, batch const& other, requires_arch) noexcept { - if (std::is_signed::value) + // Saturating sub for 8/16-bit integers needs AVX512BW; on AVX512F + // fall back to the AVX2 implementation on each 256-bit half. + XSIMD_IF_CONSTEXPR(sizeof(T) == 1 || sizeof(T) == 2) + { + return detail::fwd_to_avx([](__m256i s, __m256i o) noexcept + { return ssub(batch(s), batch(o), avx2 {}); }, + self, other); + } + else if (std::is_signed::value) { - return sadd(self, -other); + auto mask = other < 0; + auto self_overflow_branch = min(std::numeric_limits::max() + other, self); + auto self_underflow_branch = max(std::numeric_limits::min() + other, self); + return select(mask, self_overflow_branch, self_underflow_branch) - other; } else { diff --git a/include/xsimd/arch/xsimd_scalar.hpp b/include/xsimd/arch/xsimd_scalar.hpp index 28d0ba21e..044b24844 100644 --- a/include/xsimd/arch/xsimd_scalar.hpp +++ b/include/xsimd/arch/xsimd_scalar.hpp @@ -720,7 +720,21 @@ namespace xsimd { if (std::numeric_limits::is_signed) { - return sadd(lhs, (T)-rhs); + // Compute the saturating difference directly. Using sadd(lhs, -rhs) + // is wrong when rhs == numeric_limits::min(), because -rhs is not + // representable. + if ((rhs < 0) && (lhs > std::numeric_limits::max() + rhs)) + { + return std::numeric_limits::max(); + } + else if ((rhs > 0) && (lhs < std::numeric_limits::lowest() + rhs)) + { + return std::numeric_limits::lowest(); + } + else + { + return lhs - rhs; + } } else { diff --git a/test/test_xsimd_api.cpp b/test/test_xsimd_api.cpp index 454c360c9..c5e5dd8c2 100644 --- a/test/test_xsimd_api.cpp +++ b/test/test_xsimd_api.cpp @@ -500,6 +500,25 @@ TEST_CASE_TEMPLATE("[xsimd api | integral types functions]", B, INTEGRAL_TYPES) } } +// Subtracting the type minimum must saturate, not wrap. Previously this was +// computed as sadd(x, -min), and -min is not representable. +TEST_CASE_TEMPLATE("[xsimd api | ssub at type minimum]", B, INTEGRAL_TYPES) +{ + using value_type = typename scalar_type::type; + value_type lo = std::numeric_limits::min(); + value_type hi = std::numeric_limits::max(); + if (std::numeric_limits::is_signed) + { + CHECK_EQ(extract(xsimd::ssub(B(value_type(0)), B(lo))), hi); + CHECK_EQ(extract(xsimd::ssub(B(value_type(1)), B(lo))), hi); + } + else + { + // lo == 0 for unsigned types. + CHECK_EQ(extract(xsimd::ssub(B(value_type(5)), B(lo))), value_type(5)); + } +} + /* * Functions that apply on floating points types only */