From 986bfcbf8a4778ce33323a346a5f57c549731379 Mon Sep 17 00:00:00 2001 From: Monish Peddapally <93381641+monishpeddapally@users.noreply.github.com> Date: Mon, 12 Oct 2026 17:59:07 -0400 Subject: [PATCH] Fix nearbyint and rint in directed rounding modes The common kernel rounds the absolute input before restoring its sign, which reverses directed rounding for negative values. Keep the input signed during addition and subtraction, and copy its sign to the result to preserve signed zero. Use the current rounding direction in SSE4.1, AVX and AVX-512 instructions instead of selecting round-to-nearest. Add mixed-lane regressions for float and double across the four standard rounding modes. Fixes #1424 Signed-off-by: Monish Peddapally <93381641+monishpeddapally@users.noreply.github.com> --- .../xsimd/arch/common/xsimd_common_math.hpp | 9 +-- include/xsimd/arch/xsimd_avx.hpp | 4 +- include/xsimd/arch/xsimd_avx512f.hpp | 4 +- include/xsimd/arch/xsimd_sse4_1.hpp | 4 +- test/test_rounding.cpp | 58 +++++++++++++++++++ 5 files changed, 69 insertions(+), 10 deletions(-) diff --git a/include/xsimd/arch/common/xsimd_common_math.hpp b/include/xsimd/arch/common/xsimd_common_math.hpp index ad8a2741a..41a4e6dd5 100644 --- a/include/xsimd/arch/common/xsimd_common_math.hpp +++ b/include/xsimd/arch/common/xsimd_common_math.hpp @@ -1906,10 +1906,11 @@ namespace xsimd batch_type s = bitofsign(self); batch_type v = self ^ s; batch_type t2n = constants::twotonmb(); - batch_type d0 = v + t2n; - detail::reassociation_barrier(d0, "prevent collapsing (v + 2^n) - 2^n back to v"); - batch_type d = d0 - t2n; - return s ^ select(v < t2n, d, v); + batch_type signed_t2n = t2n | s; + batch_type d0 = self + signed_t2n; + detail::reassociation_barrier(d0, "prevent collapsing (x + c) - c back to x"); + batch_type d = copysign(d0 - signed_t2n, self); + return select(v < t2n, d, self); } } template diff --git a/include/xsimd/arch/xsimd_avx.hpp b/include/xsimd/arch/xsimd_avx.hpp index f06a603e4..e634aebca 100644 --- a/include/xsimd/arch/xsimd_avx.hpp +++ b/include/xsimd/arch/xsimd_avx.hpp @@ -1345,12 +1345,12 @@ namespace xsimd template XSIMD_INLINE batch nearbyint(batch const& self, requires_arch) noexcept { - return _mm256_round_ps(self, _MM_FROUND_TO_NEAREST_INT); + return _mm256_round_ps(self, _MM_FROUND_CUR_DIRECTION); } template XSIMD_INLINE batch nearbyint(batch const& self, requires_arch) noexcept { - return _mm256_round_pd(self, _MM_FROUND_TO_NEAREST_INT); + return _mm256_round_pd(self, _MM_FROUND_CUR_DIRECTION); } // nearbyint_as_int diff --git a/include/xsimd/arch/xsimd_avx512f.hpp b/include/xsimd/arch/xsimd_avx512f.hpp index 06da6218e..0e2915d67 100644 --- a/include/xsimd/arch/xsimd_avx512f.hpp +++ b/include/xsimd/arch/xsimd_avx512f.hpp @@ -1868,12 +1868,12 @@ namespace xsimd template XSIMD_INLINE batch nearbyint(batch const& self, requires_arch) noexcept { - return _mm512_roundscale_round_ps(self, _MM_FROUND_TO_NEAREST_INT, _MM_FROUND_CUR_DIRECTION); + return _mm512_roundscale_round_ps(self, _MM_FROUND_CUR_DIRECTION, _MM_FROUND_CUR_DIRECTION); } template XSIMD_INLINE batch nearbyint(batch const& self, requires_arch) noexcept { - return _mm512_roundscale_round_pd(self, _MM_FROUND_TO_NEAREST_INT, _MM_FROUND_CUR_DIRECTION); + return _mm512_roundscale_round_pd(self, _MM_FROUND_CUR_DIRECTION, _MM_FROUND_CUR_DIRECTION); } // nearbyint_as_int diff --git a/include/xsimd/arch/xsimd_sse4_1.hpp b/include/xsimd/arch/xsimd_sse4_1.hpp index 73af32cd6..27fe585c3 100644 --- a/include/xsimd/arch/xsimd_sse4_1.hpp +++ b/include/xsimd/arch/xsimd_sse4_1.hpp @@ -408,12 +408,12 @@ namespace xsimd template XSIMD_INLINE batch nearbyint(batch const& self, requires_arch) noexcept { - return _mm_round_ps(self, _MM_FROUND_TO_NEAREST_INT); + return _mm_round_ps(self, _MM_FROUND_CUR_DIRECTION); } template XSIMD_INLINE batch nearbyint(batch const& self, requires_arch) noexcept { - return _mm_round_pd(self, _MM_FROUND_TO_NEAREST_INT); + return _mm_round_pd(self, _MM_FROUND_CUR_DIRECTION); } // select diff --git a/test/test_rounding.cpp b/test/test_rounding.cpp index 681aaaef7..d6177faf8 100644 --- a/test/test_rounding.cpp +++ b/test/test_rounding.cpp @@ -14,6 +14,8 @@ #include "test_utils.hpp" +#include + template struct rounding_test { @@ -167,4 +169,60 @@ TEST_CASE_TEMPLATE("[rounding]", B, BATCH_FLOAT_TYPES) rounding_test Test; Test.test_rounding_functions(); } + +#ifndef __FAST_MATH__ +TEST_CASE_TEMPLATE("[rounding] current rounding mode", B, BATCH_FLOAT_TYPES) +{ + using value_type = typename B::value_type; + struct rounding_mode_guard + { + int saved = std::fegetround(); + ~rounding_mode_guard() { std::fesetround(saved); } + } guard; + + const std::array input = { + 0., -0., 0.25, -0.25, 0.5, -0.5, 1.25, -1.25, + 1.75, -1.75, 2.5, -2.5, 3.5, -3.5, 2., -2. + }; + const std::array modes = { FE_TONEAREST, FE_DOWNWARD, FE_UPWARD, FE_TOWARDZERO }; + const value_type expected[][16] = { + { 0., -0., 0., -0., 0., -0., 1., -1., 2., -2., 2., -2., 4., -4., 2., -2. }, + { 0., -0., 0., -1., 0., -1., 1., -2., 1., -2., 2., -3., 3., -4., 2., -2. }, + { 0., -0., 1., -0., 1., -0., 2., -1., 2., -1., 3., -2., 4., -3., 2., -2. }, + { 0., -0., 0., -0., 0., -0., 1., -1., 1., -1., 2., -2., 3., -3., 2., -2. } + }; + + for (size_t mode = 0; mode < modes.size(); ++mode) + { + if (std::fesetround(modes[mode]) != 0) + { + continue; + } + REQUIRE(std::fegetround() == modes[mode]); + for (size_t offset = 0; offset < input.size(); offset += B::size) + { + std::array values, nearby, rounded; + for (size_t lane = 0; lane < B::size; ++lane) + { + // Reload after changing the mode so arithmetic cannot be hoisted. + volatile value_type value = input[(offset + lane) % input.size()]; + values[lane] = value; + } + B batch = B::load_unaligned(values.data()); + xsimd::nearbyint(batch).store_unaligned(nearby.data()); + xsimd::rint(batch).store_unaligned(rounded.data()); + for (size_t lane = 0; lane < B::size; ++lane) + { + CAPTURE(modes[mode]); + CAPTURE(values[lane]); + value_type reference = expected[mode][(offset + lane) % input.size()]; + CHECK(nearby[lane] == reference); + CHECK(rounded[lane] == reference); + CHECK(std::signbit(nearby[lane]) == std::signbit(reference)); + CHECK(std::signbit(rounded[lane]) == std::signbit(reference)); + } + } + } +} +#endif #endif