From 4bc84b22f94ce55a1a6462f858f87e5341ccf7a9 Mon Sep 17 00:00:00 2001 From: Demetrios Agourakis Date: Fri, 30 Oct 2026 10:41:39 -0300 Subject: [PATCH 1/2] Handle signed zero pairs without invalid division in atan2 --- .../xsimd/arch/common/xsimd_common_trigo.hpp | 11 ++++- test/test_complex_exponential.cpp | 34 +++++++++++++++ test/test_trigonometric.cpp | 41 +++++++++++++++++++ 3 files changed, 84 insertions(+), 2 deletions(-) diff --git a/include/xsimd/arch/common/xsimd_common_trigo.hpp b/include/xsimd/arch/common/xsimd_common_trigo.hpp index 00ad75076..ad8bfaf7c 100644 --- a/include/xsimd/arch/common/xsimd_common_trigo.hpp +++ b/include/xsimd/arch/common/xsimd_common_trigo.hpp @@ -358,11 +358,18 @@ namespace xsimd const batch_type q = abs(self / other); const batch_type q_p = abs(other / self); #else - const batch_type q = abs(self / other); + const auto both_zero = (self == batch_type(0.)) && (other == batch_type(0.)); + const batch_type q = abs(select(both_zero, batch_type(1.), self) + / select(both_zero, batch_type(1.), other)); const batch_type q_p = 1. / q; #endif const batch_type z = detail::kernel_atan(q, q_p); - return select(other > batch_type(0.), z, constants::pi() - z) * signnz(self); + batch_type result = select(other > batch_type(0.), z, constants::pi() - z) * signnz(self); +#ifndef __FAST_MATH__ + const auto zero_angle = select(signnz(other) < batch_type(0.), constants::pi(), batch_type(0.)); + result = select(both_zero, copysign(zero_angle, self), result); +#endif + return result; } // cos diff --git a/test/test_complex_exponential.cpp b/test/test_complex_exponential.cpp index 8dffdf5c4..bb128911c 100644 --- a/test/test_complex_exponential.cpp +++ b/test/test_complex_exponential.cpp @@ -213,4 +213,38 @@ TEST_CASE_TEMPLATE("[complex exponential]", B, BATCH_COMPLEX_TYPES) Test.test_sign(); } } +#ifndef __FAST_MATH__ +TEST_CASE_TEMPLATE("[complex log zero angle]", 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 samples = { + value_type(zero, zero), value_type(zero, -zero), + value_type(-zero, zero), value_type(-zero, -zero) + }; + std::array inputs; + std::array results; + for (size_t start = 0; start < samples.size(); ++start) + { + for (size_t lane = 0; lane < B::size; ++lane) + { + inputs[lane] = samples[(start + lane) % samples.size()]; + } + xsimd::log(B::load_unaligned(inputs.data())).store_unaligned(results.data()); + for (size_t lane = 0; lane < B::size; ++lane) + { + const auto expected = std::log(inputs[lane]); + CHECK(std::isinf(results[lane].real())); + CHECK(std::signbit(results[lane].real())); + CHECK_EQ(results[lane].imag(), doctest::Approx(expected.imag()).epsilon(real_type(4) * std::numeric_limits::epsilon())); + CHECK_EQ(std::signbit(results[lane].imag()), std::signbit(expected.imag())); + if (expected.imag() == zero) + { + CHECK_EQ(results[lane].imag(), zero); + } + } + } +} +#endif #endif diff --git a/test/test_trigonometric.cpp b/test/test_trigonometric.cpp index d3e8862a6..364c7861d 100644 --- a/test/test_trigonometric.cpp +++ b/test/test_trigonometric.cpp @@ -194,4 +194,45 @@ TEST_CASE_TEMPLATE("[trigonometric]", B, BATCH_FLOAT_TYPES) Test.test_reciprocal_functions(); } } +#ifndef __FAST_MATH__ +TEST_CASE_TEMPLATE("[atan2 signed zero pairs]", B, BATCH_FLOAT_TYPES) +{ + using value_type = typename B::value_type; + const value_type zero(0); + const std::array, 20> samples = { + std::pair { zero, zero }, std::pair { -zero, zero }, + std::pair { zero, -zero }, std::pair { -zero, -zero }, + std::pair { zero, value_type(1) }, std::pair { -zero, value_type(1) }, + std::pair { zero, value_type(-1) }, std::pair { -zero, value_type(-1) }, + std::pair { value_type(1), zero }, std::pair { value_type(-1), zero }, + std::pair { value_type(1), -zero }, std::pair { value_type(-1), -zero }, + std::pair { value_type(1), value_type(1) }, std::pair { value_type(-1), value_type(1) }, + std::pair { value_type(1), value_type(-1) }, std::pair { value_type(-1), value_type(-1) }, + std::pair { value_type(3), value_type(4) }, std::pair { value_type(-3), value_type(4) }, + std::pair { value_type(3), value_type(-4) }, std::pair { value_type(-3), value_type(-4) } + }; + std::array ys, xs, results; + const auto tolerance = value_type(64) * std::numeric_limits::epsilon(); + for (size_t start = 0; start < samples.size(); ++start) + { + for (size_t lane = 0; lane < B::size; ++lane) + { + const auto sample = samples[(start + lane) % samples.size()]; + ys[lane] = sample.first; + xs[lane] = sample.second; + } + xsimd::atan2(B::load_unaligned(ys.data()), B::load_unaligned(xs.data())).store_unaligned(results.data()); + for (size_t lane = 0; lane < B::size; ++lane) + { + const auto expected = std::atan2(ys[lane], xs[lane]); + CHECK_EQ(results[lane], doctest::Approx(expected).epsilon(tolerance).scale(1)); + CHECK_EQ(std::signbit(results[lane]), std::signbit(expected)); + if (expected == zero) + { + CHECK_EQ(results[lane], expected); + } + } + } +} +#endif #endif From 94a8f71a3a57e194a9f3d18e2dedfbfc74955dca Mon Sep 17 00:00:00 2001 From: Demetrios Agourakis Date: Fri, 30 Oct 2026 12:44:49 -0300 Subject: [PATCH 2/2] Handle pairs of infinite arguments in atan2 --- .../xsimd/arch/common/xsimd_common_trigo.hpp | 8 ++- test/test_trigonometric.cpp | 49 +++++++++++++++++++ 2 files changed, 55 insertions(+), 2 deletions(-) diff --git a/include/xsimd/arch/common/xsimd_common_trigo.hpp b/include/xsimd/arch/common/xsimd_common_trigo.hpp index ad8bfaf7c..462c94026 100644 --- a/include/xsimd/arch/common/xsimd_common_trigo.hpp +++ b/include/xsimd/arch/common/xsimd_common_trigo.hpp @@ -359,8 +359,12 @@ namespace xsimd const batch_type q_p = abs(other / self); #else const auto both_zero = (self == batch_type(0.)) && (other == batch_type(0.)); - const batch_type q = abs(select(both_zero, batch_type(1.), self) - / select(both_zero, batch_type(1.), other)); + auto use_unit_ratio = both_zero; +#ifndef XSIMD_NO_INFINITIES + use_unit_ratio = use_unit_ratio || (isinf(self) && isinf(other)); +#endif + const batch_type q = abs(select(use_unit_ratio, batch_type(1.), self) + / select(use_unit_ratio, batch_type(1.), other)); const batch_type q_p = 1. / q; #endif const batch_type z = detail::kernel_atan(q, q_p); diff --git a/test/test_trigonometric.cpp b/test/test_trigonometric.cpp index 364c7861d..bb91cbfd1 100644 --- a/test/test_trigonometric.cpp +++ b/test/test_trigonometric.cpp @@ -234,5 +234,54 @@ TEST_CASE_TEMPLATE("[atan2 signed zero pairs]", B, BATCH_FLOAT_TYPES) } } } +#if !defined(XSIMD_NO_INFINITIES) && !defined(XSIMD_NO_NANS) +TEST_CASE_TEMPLATE("[atan2 infinite argument pairs]", B, BATCH_FLOAT_TYPES) +{ + using value_type = typename B::value_type; + const value_type zero(0); + const auto inf = std::numeric_limits::infinity(); + const auto nan = std::numeric_limits::quiet_NaN(); + const std::array, 20> samples = { + std::pair { inf, inf }, std::pair { -inf, inf }, + std::pair { inf, -inf }, std::pair { -inf, -inf }, + std::pair { inf, value_type(1) }, std::pair { -inf, value_type(1) }, + std::pair { inf, value_type(-1) }, std::pair { -inf, value_type(-1) }, + std::pair { value_type(1), inf }, std::pair { value_type(-1), inf }, + std::pair { value_type(1), -inf }, std::pair { value_type(-1), -inf }, + std::pair { zero, inf }, std::pair { -zero, inf }, + std::pair { zero, -inf }, std::pair { -zero, -inf }, + std::pair { nan, value_type(1) }, std::pair { value_type(1), nan }, + std::pair { value_type(1), value_type(1) }, std::pair { value_type(-1), value_type(-1) } + }; + std::array ys, xs, results; + const auto tolerance = value_type(64) * std::numeric_limits::epsilon(); + for (size_t start = 0; start < samples.size(); ++start) + { + for (size_t lane = 0; lane < B::size; ++lane) + { + const auto sample = samples[(start + lane) % samples.size()]; + ys[lane] = sample.first; + xs[lane] = sample.second; + } + xsimd::atan2(B::load_unaligned(ys.data()), B::load_unaligned(xs.data())).store_unaligned(results.data()); + for (size_t lane = 0; lane < B::size; ++lane) + { + const auto expected = std::atan2(ys[lane], xs[lane]); + if (std::isnan(expected)) + { + CHECK(std::isnan(results[lane])); + continue; + } + CHECK(std::isfinite(results[lane])); + CHECK_EQ(results[lane], doctest::Approx(expected).epsilon(tolerance).scale(1)); + CHECK_EQ(std::signbit(results[lane]), std::signbit(expected)); + if (expected == zero) + { + CHECK_EQ(results[lane], zero); + } + } + } +} +#endif #endif #endif