Fix nearbyint and rint in directed rounding modes

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>();
-                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 <class A>

include/xsimd/arch/xsimd_avx.hpp

@@ -1345,12 +1345,12 @@ namespace xsimd
         template <class A>
         XSIMD_INLINE batch<float, A> nearbyint(batch<float, A> const& self, requires_arch<avx>) noexcept
         {
-            return _mm256_round_ps(self, _MM_FROUND_TO_NEAREST_INT);
+            return _mm256_round_ps(self, _MM_FROUND_CUR_DIRECTION);
         }
         template <class A>
         XSIMD_INLINE batch<double, A> nearbyint(batch<double, A> const& self, requires_arch<avx>) noexcept
         {
-            return _mm256_round_pd(self, _MM_FROUND_TO_NEAREST_INT);
+            return _mm256_round_pd(self, _MM_FROUND_CUR_DIRECTION);
         }
 
         // nearbyint_as_int

include/xsimd/arch/xsimd_avx512f.hpp

@@ -1868,12 +1868,12 @@ namespace xsimd
         template <class A>
         XSIMD_INLINE batch<float, A> nearbyint(batch<float, A> const& self, requires_arch<avx512f>) 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 <class A>
         XSIMD_INLINE batch<double, A> nearbyint(batch<double, A> const& self, requires_arch<avx512f>) 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

include/xsimd/arch/xsimd_sse4_1.hpp

@@ -408,12 +408,12 @@ namespace xsimd
         template <class A>
         XSIMD_INLINE batch<float, A> nearbyint(batch<float, A> const& self, requires_arch<sse4_1>) noexcept
         {
-            return _mm_round_ps(self, _MM_FROUND_TO_NEAREST_INT);
+            return _mm_round_ps(self, _MM_FROUND_CUR_DIRECTION);
         }
         template <class A>
         XSIMD_INLINE batch<double, A> nearbyint(batch<double, A> const& self, requires_arch<sse4_1>) noexcept
         {
-            return _mm_round_pd(self, _MM_FROUND_TO_NEAREST_INT);
+            return _mm_round_pd(self, _MM_FROUND_CUR_DIRECTION);
         }
 
         // select

test/test_rounding.cpp

@@ -14,6 +14,8 @@
 
 #include "test_utils.hpp"
 
+#include <cfenv>
+
 template <class B>
 struct rounding_test
 {
@@ -167,4 +169,60 @@ TEST_CASE_TEMPLATE("[rounding]", B, BATCH_FLOAT_TYPES)
     rounding_test<B> 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<value_type, 16> 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<int, 4> 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<value_type, B::size> 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