From f961a6e81462c67608e418beb557f511a4a017e1 Mon Sep 17 00:00:00 2001 From: Demetrios Chiuratto Agourakis Date: Mon, 26 Oct 2026 16:05:18 +0000 Subject: [PATCH 1/3] Preserve signed imaginary zero in complex sqrt --- .../xsimd/arch/common/xsimd_common_math.hpp | 6 +-- test/test_complex_power.cpp | 37 +++++++++++++++++++ 2 files changed, 40 insertions(+), 3 deletions(-) diff --git a/include/xsimd/arch/common/xsimd_common_math.hpp b/include/xsimd/arch/common/xsimd_common_math.hpp index ad8a2741a..8ab7e0e59 100644 --- a/include/xsimd/arch/common/xsimd_common_math.hpp +++ b/include/xsimd/arch/common/xsimd_common_math.hpp @@ -2323,9 +2323,9 @@ namespace xsimd real_batch ze(0.); return select(y == ze, - select(x == ze, - batch_type(ze, ze), - select(x < ze, batch_type(ze, sqrt_x), batch_type(sqrt_x, ze))), + select(x < ze, + batch_type(ze, copysign(sqrt_x, y)), + batch_type(sqrt_x, y)), select(x == ze, select(y > ze, batch_type(sqrt_hy, sqrt_hy), batch_type(sqrt_hy, -sqrt_hy)), resg)); diff --git a/test/test_complex_power.cpp b/test/test_complex_power.cpp index 99bb82819..cbb7dd093 100644 --- a/test/test_complex_power.cpp +++ b/test/test_complex_power.cpp @@ -256,4 +256,41 @@ TEST_CASE_TEMPLATE("[complex power]", B, BATCH_COMPLEX_TYPES) Test.test_sqrt_pp(); } } +#ifndef __FAST_MATH__ +TEST_CASE_TEMPLATE("[complex sqrt signed zero]", B, BATCH_COMPLEX_TYPES) +{ + using value_type = typename B::value_type; + using real_type = typename value_type::value_type; + const real_type zero(0); + const std::array real_parts = { + -std::numeric_limits::max(), real_type(-4), real_type(-1), + -std::numeric_limits::denorm_min(), -zero, zero, + std::numeric_limits::denorm_min(), real_type(1), real_type(4), + std::numeric_limits::max() + }; + std::array inputs; + std::array results; + + for (size_t start = 0; start < 2 * real_parts.size(); ++start) + { + for (size_t lane = 0; lane < B::size; ++lane) + { + const size_t index = (start + lane) % (2 * real_parts.size()); + inputs[lane] = value_type(real_parts[index / 2], index % 2 ? -zero : zero); + } + + xsimd::sqrt(B::load_unaligned(inputs.data())).store_unaligned(results.data()); + + for (size_t lane = 0; lane < B::size; ++lane) + { + const auto expected = std::sqrt(inputs[lane]); + CHECK_EQ(results[lane].real(), expected.real()); + CHECK_EQ(results[lane].imag(), expected.imag()); + CHECK_EQ(std::signbit(results[lane].real()), std::signbit(expected.real())); + CHECK_EQ(std::signbit(results[lane].imag()), std::signbit(expected.imag())); + } + } +} +#endif + #endif From 72d4e6f3e2e6965e5bc1d72258ead17b45a5d97f Mon Sep 17 00:00:00 2001 From: Demetrios Agourakis Date: Wed, 28 Oct 2026 06:18:52 -0300 Subject: [PATCH 2/3] Use real-axis identity in complex sqrt signed-zero regression --- test/test_complex_power.cpp | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/test/test_complex_power.cpp b/test/test_complex_power.cpp index cbb7dd093..557bef1d3 100644 --- a/test/test_complex_power.cpp +++ b/test/test_complex_power.cpp @@ -283,7 +283,15 @@ TEST_CASE_TEMPLATE("[complex sqrt signed zero]", B, BATCH_COMPLEX_TYPES) for (size_t lane = 0; lane < B::size; ++lane) { - const auto expected = std::sqrt(inputs[lane]); + // Use the real-axis identity: some complex sqrt implementations + // evaluate polar coordinates and leave a nonzero real residual + // on the negative real axis. + const auto x = inputs[lane].real(); + const auto y = inputs[lane].imag(); + const auto magnitude = std::sqrt(std::abs(x)); + const value_type expected = x < zero + ? value_type(zero, std::copysign(magnitude, y)) + : value_type(magnitude, y); CHECK_EQ(results[lane].real(), expected.real()); CHECK_EQ(results[lane].imag(), expected.imag()); CHECK_EQ(std::signbit(results[lane].real()), std::signbit(expected.real())); From acdd9f4ea9d16fa1abbf9bdf94949d5025fd58e1 Mon Sep 17 00:00:00 2001 From: Demetrios Agourakis Date: Thu, 29 Oct 2026 06:18:31 -0300 Subject: [PATCH 3/3] Fix complex sqrt signed-zero test oracle portability Use the real SIMD backend for magnitude expectations while keeping exact component and sign checks. This respects backend rounding and underflow behavior without changing the complex sqrt implementation. --- test/test_complex_power.cpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/test/test_complex_power.cpp b/test/test_complex_power.cpp index 557bef1d3..0e00b4327 100644 --- a/test/test_complex_power.cpp +++ b/test/test_complex_power.cpp @@ -270,6 +270,7 @@ TEST_CASE_TEMPLATE("[complex sqrt signed zero]", B, BATCH_COMPLEX_TYPES) }; std::array inputs; std::array results; + std::array magnitudes; for (size_t start = 0; start < 2 * real_parts.size(); ++start) { @@ -279,16 +280,17 @@ TEST_CASE_TEMPLATE("[complex sqrt signed zero]", B, BATCH_COMPLEX_TYPES) inputs[lane] = value_type(real_parts[index / 2], index % 2 ? -zero : zero); } - xsimd::sqrt(B::load_unaligned(inputs.data())).store_unaligned(results.data()); + const auto batch_inputs = B::load_unaligned(inputs.data()); + xsimd::sqrt(batch_inputs).store_unaligned(results.data()); + // Test complex composition while preserving the real backend's + // rounding and underflow behavior. + xsimd::sqrt(xsimd::abs(batch_inputs.real())).store_unaligned(magnitudes.data()); for (size_t lane = 0; lane < B::size; ++lane) { - // Use the real-axis identity: some complex sqrt implementations - // evaluate polar coordinates and leave a nonzero real residual - // on the negative real axis. const auto x = inputs[lane].real(); const auto y = inputs[lane].imag(); - const auto magnitude = std::sqrt(std::abs(x)); + const auto magnitude = magnitudes[lane]; const value_type expected = x < zero ? value_type(zero, std::copysign(magnitude, y)) : value_type(magnitude, y);