Skip to content

Preserve signed imaginary zero in complex sqrt #1433

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
agourakis82 wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Preserve signed imaginary zero in complex sqrt #1433

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
f961a6e
Preserve signed imaginary zero in complex sqrt
agourakis82 Oct 26, 2026
72d4e6f
Use real-axis identity in complex sqrt signed-zero regression
agourakis82 Oct 28, 2026
acdd9f4
Fix complex sqrt signed-zero test oracle portability
agourakis82 Oct 29, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions include/xsimd/arch/common/xsimd_common_math.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,

@serge-sans-paille serge-sans-paille Oct 28, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I like how this also should make the code denser by avoiding a select.
The logic looks good too.

batch_type(ze, copysign(sqrt_x, y)),

@serge-sans-paille serge-sans-paille Oct 28, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

mmmh, that part looks strange, why do you change the sign of sqrt_x here? I don't see why sqrt_x would be zero...

@agourakis82 agourakis82 Oct 29, 2026

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here sqrt_x is sqrt(abs(x)). In the y == 0 && x < 0 branch, it becomes the magnitude of the imaginary component, so it does not need to be zero.

For example, sqrt(complex(-4, -0)) should return (0, -2), whereas sqrt(complex(-4, +0)) should return (0, +2). The copysign(sqrt_x, y) selects the corresponding side of the negative-real-axis branch cut. Without it, both cases would return a positive imaginary component.

The change preserves the nonnegative real component; it applies the sign of y to the imaginary component in this branch.

batch_type(sqrt_x, y)),
select(x == ze,
select(y > ze, batch_type(sqrt_hy, sqrt_hy), batch_type(sqrt_hy, -sqrt_hy)),
resg));
Expand Down
47 changes: 47 additions & 0 deletions test/test_complex_power.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -256,4 +256,51 @@ 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_type, 10> real_parts = {
-std::numeric_limits<real_type>::max(), real_type(-4), real_type(-1),
-std::numeric_limits<real_type>::denorm_min(), -zero, zero,
std::numeric_limits<real_type>::denorm_min(), real_type(1), real_type(4),
std::numeric_limits<real_type>::max()
};
std::array<value_type, B::size> inputs;
std::array<value_type, B::size> results;
std::array<real_type, B::size> magnitudes;

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);
}

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)
{
const auto x = inputs[lane].real();
const auto y = inputs[lane].imag();
const auto magnitude = magnitudes[lane];
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()));
CHECK_EQ(std::signbit(results[lane].imag()), std::signbit(expected.imag()));
}
}
}
#endif

#endif
Loading