Skip to content

Preserve signed imaginary zero in complex sqrt - #1433

Open
agourakis82 wants to merge 3 commits into
xtensor-stack:masterfrom
agourakis82:fix/complex-sqrt-signed-zero-1432
Open

agourakis82 wants to merge 3 commits into
xtensor-stack:masterfrom
agourakis82:fix/complex-sqrt-signed-zero-1432

Conversation

@agourakis82

@agourakis82 agourakis82 commented Oct 26, 2026 •

Copy link
Copy Markdown

For finite real inputs on the real axis, xsimd::sqrt discards the sign of the imaginary zero. For example, sqrt(complex(-4, -0)) returns (0, 2) instead of (0, -2), selecting the wrong side of the branch cut. Positive real inputs and complex zero also lose an imaginary -0.

A standalone consumer also demonstrates a value-level consequence: composing the vectorized root with scalar std::exp(i * root) for (-4, -0) produces approximately 0.135335 before the fix, whereas the scalar reference produces 7.389056. After the fix the composition agrees with the reference. This is a synthetic composition that isolates the sqrt defect, not a claim of an observed failure in a downstream application.

Preserve the imaginary zero for nonnegative real inputs and use copysign for the imaginary root of negative real inputs. The real component remains nonnegative, including for a negative real zero.

The regression checks both components and their sign bits using the real-axis identities to construct the expected result, rather than complex std::sqrt. It uses the real SIMD backend for the magnitude and independently constructs component placement and zero signs. This checks complex composition while respecting the real backend's rounding and underflow behavior; real sqrt accuracy is covered separately. Mixed SIMD lanes cover positive and negative normal values, the smallest subnormal, the largest finite value, and both real zeros, with both signs of imaginary zero. The test is disabled under FAST_MATH, where signed-zero behavior is not guaranteed.

This addresses the finite real-axis sqrt cases from #1432. The other functions and non-real-axis cases reported there remain outside this change; this PR should not close the issue.

Native validation with GCC 13.3, C++23 and doctest 2.4.9:

  • The new regression fails against the unchanged base for both float and double in SSE2 and AVX2 builds, then passes with the fix.
  • All tests in test_complex_power.cpp pass in SSE2 and AVX2 builds.
  • The same unit passes AddressSanitizer and UndefinedBehaviorSanitizer.
  • The full CMake/CTest suite passes with TARGET_ARCH=sse2: 434 doctest cases, 78,104 assertions; CTest 1/1.
  • The original issue reproducer was run unchanged before and after the patch in both SIMD configurations; the two sqrt examples now agree with std::sqrt, while the other reported failures remain visible.

Portability follow-up (C++23, AppleClang, native ARM64): the updated regression passes 2 cases / 480 assertions; the complete complex-power unit passes 4 cases / 486 assertions. Independent runs with flush-to-zero off and on pass. Complete private-header controls removing each sign correction are rejected, with no tolerance increase or input skip. These are bounded native checks, not an ARMv7 or i386 replay.

These are native Linux x86-64 results. Other architectures and official CI remain to be checked.

AI assistance was used to investigate the report, prepare the patch, and run validation.

select(x == ze,
batch_type(ze, ze),
select(x < ze, batch_type(ze, sqrt_x), batch_type(sqrt_x, ze))),
select(x < ze,

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.

@serge-sans-paille serge-sans-paille left a comment

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.

CI failing ^^!

@agourakis82

Copy link
Copy Markdown
Author

Thanks for pointing this out. In 72d4e6f, I changed the regression reference to use the real-axis identities rather than complex std::sqrt, and I've updated the PR description to match.

The workflows for this commit are currently waiting for approval. Could you approve them so we can check the updated test across the CI configurations?

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)),

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...

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.

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.
@agourakis82

Copy link
Copy Markdown
Author

The updated run has 86 passing checks and two failures in the MSYS2 jobs. MINGW32 stopped during environment setup because the package mirror timed out. CLANG64 reached the test invocation but exited with code 127, without reporting test results; I haven't established the cause of that exit yet. Could you rerun those two jobs?

@agourakis82

Copy link
Copy Markdown
Author

I captured the native Windows status in a fresh paired CLANG64 run: both master ebe9991 and this PR's merge d46d757 start successfully, but their full suites raise STATUS_FLOAT_INVALID_OPERATION (0xC0000090), which MSYS2 presents as 127. This is not a DLL-loading failure in that reproduction.

The failure localizes to the existing [power] tests. A diagnostic copy records the fault in CHECK_BATCH_EQ; disassembly identifies a speculative vdivps inside the shared safe_division comparison helper. The power tests unmask invalid exceptions, while the helper was compiled assuming exceptions remain masked. The focused sqrt regression passes independently: 2 cases / 960 assertions.

I opened #1440 with a small shared test-harness repair: scoped Clang exception semantics, preserving the guards, arithmetic, assertions and invalid-exception unmasking. No sqrt implementation changes are included there.

Fresh full C++23 Windows qualification reproduces the original failures before the repair and passes afterward: master 430 cases / 151245 assertions; this PR's merge 432 / 152205. The full runs exclude no tests, and the original MSYS launch returns zero. These fork results are separate from official CI, and the historical job's original native status remains unavailable. This PR's published source is unchanged by the diagnosis.

Prepared with OpenAI Codex assistance.

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants