From 6513d415150831aebba655c8b5919cc1e51bc9de Mon Sep 17 00:00:00 2001 From: Fikret Ardal Date: Wed, 11 Nov 2026 16:51:07 -0500 Subject: [PATCH 1/2] Fix complex asin, acos, tan, tanh and log1p mismatches with std Fix #1432: asin/acos/asinh/acosh use Kahan's formulas, tan/tanh use FreeBSD's ctanh, log1p keeps the sign of a -0 imaginary part, and atan2(+-0, +-0) no longer returns NaN. Add tests with every input from the issue. Assisted-By: Claude Opus 5.5 --- .../xsimd/arch/common/xsimd_common_math.hpp | 12 +- .../xsimd/arch/common/xsimd_common_trigo.hpp | 146 +++++++++++++----- include/xsimd/arch/xsimd_scalar.hpp | 4 +- test/test_complex_exponential.cpp | 11 ++ test/test_complex_trigonometric.cpp | 36 +++++ test/test_trigonometric.cpp | 13 ++ test/test_utils.hpp | 34 ++++ 7 files changed, 210 insertions(+), 46 deletions(-) diff --git a/include/xsimd/arch/common/xsimd_common_math.hpp b/include/xsimd/arch/common/xsimd_common_math.hpp index 8ab7e0e59..b9484fa81 100644 --- a/include/xsimd/arch/common/xsimd_common_math.hpp +++ b/include/xsimd/arch/common/xsimd_common_math.hpp @@ -1875,11 +1875,13 @@ namespace xsimd using real_batch = typename batch_type::real_batch; batch_type u = 1 + self; batch_type logu = log(u); - return select(u == batch_type(1.), - self, - select(u.real() <= real_batch(0.), - logu, - logu * self / (u - batch_type(1.)))); + batch_type res = select(u == batch_type(1.), + self, + select(u.real() <= real_batch(0.), + logu, + logu * self / (u - batch_type(1.)))); + // Im(log(1 + z)) = arg(1 + z) has the sign of Im(z). Restore it. + return { res.real(), copysign(res.imag(), self.imag()) }; } // mod diff --git a/include/xsimd/arch/common/xsimd_common_trigo.hpp b/include/xsimd/arch/common/xsimd_common_trigo.hpp index 00ad75076..6119c65fb 100644 --- a/include/xsimd/arch/common/xsimd_common_trigo.hpp +++ b/include/xsimd/arch/common/xsimd_common_trigo.hpp @@ -46,13 +46,70 @@ namespace xsimd x = select(self < batch_type(-0.5), constants::pi() - x, x); return select(x_larger_05, x, constants::pio2() - x); } + namespace detail + { + // For z = x + iy, returns Re sqrt(1 - z^2) as the real part and |Im asin(z)| as the imaginary part. + // With p = sqrt(1 + z) and q = sqrt(1 - z): + // Re sqrt(1 - z^2) = Re(q * p), |Im asin(z)| = asinh(Im(conj(q) p)) + // as in Julia's asin(z::Complex): https://github.com/JuliaLang/julia/blob/master/base/complex.jl + template + XSIMD_INLINE batch, A> asin_kernel(const batch, A>& z) noexcept + { + using batch_type = batch, A>; + using real_batch = typename batch_type::real_batch; + const real_batch one(1.); + const real_batch half(0.5); + const real_batch x = abs(z.real()); + const real_batch y = abs(z.imag()); + const real_batch hy = half * y; // y / 2 + const real_batch xm1 = x - one; // x - 1 + const real_batch xp1 = x + one; // x + 1 + + // p = sqrt(1 + z) = a + bi, and q = sqrt(1 - z) with a bigger part c and a smaller part d: + // a = sqrt(((1 + x) + sqrt((1 + x)^2 + y^2)) / 2), b = y / (2a) + // c = sqrt((sqrt((1 - x)^2 + y^2) + |1 - x|) / 2), d = y / (2c) + // 1 + x >= 1, so nothing cancels in a and a >= 1. Using |1 - x| instead of 1 - x adds two + // non-negative terms, so nothing cancels in c. For x <= 1, q = c - di (c is the real part); for x > 1, + // c and d swap places: q = d - ci. + // s = sqrt((1 - x)^2 + y^2). At x = 1 it is y, taken directly: hypot squares y, and y^2 underflows + // for tiny y. + const real_batch s = select(xm1 == real_batch(0.), y, hypot(xm1, y)); + const real_batch a = sqrt(half * (hypot(xp1, y) + xp1)); + const real_batch b = hy / a; + const real_batch c = sqrt(half * (s + abs(xm1))); + const real_batch d = hy / max(c, constants::smallestposval()); // c is 0 at z = +-1 + const auto inside = x <= one; + + // Re sqrt(1 - z^2) = Re(q * p). + // x <= 1: (c - di)(a + bi) has real part ac + bd. x > 1: (d - ci)(a + bi) has real part ad + bc. + real_batch re = select(inside, fma(a, c, b * d), fma(a, d, b * c)); + + // |Im asin(z)| = asinh(Im(conj(q) p)). + // x <= 1: (c + di)(a + bi) has imaginary part ad + bc. x > 1: (d + ci)(a + bi) has imaginary part ac + bd. + const real_batch sh = select(inside, fma(a, d, b * c), fma(a, c, b * d)); + real_batch v = asinh(sh); + + // For huge x or y, (1 +- x)^2 + y^2 overflows. There Re sqrt(1 - z^2) = y and + // |Im asin(z)| = log(2 sqrt(x^2 + y^2)) up to rounding, as in FreeBSD's casinh: + // https://github.com/freebsd/freebsd-src/blob/main/lib/msun/src/catrig.c + const auto y_larger = x < y; + const real_batch m = select(y_larger, y, x); // max(|x|, |y|) + const auto huge = m >= constants::oneotwoeps(); + if (any(huge)) + { + const real_batch r = select(y_larger, x / y, y / x); // min(|x|, |y|) / max(|x|, |y|) + re = select(huge, y, re); // Re sqrt(1 - z^2) = |y| + v = select(huge, log(m) + constants::log_2() + half * log1p(r * r), v); // log(2 sqrt(x^2 + y^2)) = log(m) + log(2) + log(1 + r^2) / 2 + } + return { re, v }; + } + } template XSIMD_INLINE batch, A> acos(const batch, A>& z, requires_arch) noexcept { using batch_type = batch, A>; - using real_batch = typename batch_type::real_batch; - batch_type tmp = asin(z); - return { constants::pio2() - tmp.real(), -tmp.imag() }; + const batch_type w = detail::asin_kernel(z); // (Re sqrt(1 - z^2), |Im asin(z)|) + return { atan2(w.real(), z.real()), copysign(w.imag(), -z.imag()) }; } // acosh @@ -79,9 +136,8 @@ namespace xsimd XSIMD_INLINE batch, A> acosh(const batch, A>& z, requires_arch) noexcept { using batch_type = batch, A>; - batch_type w = acos(z); - w = batch_type(-w.imag(), w.real()); - return w; + const batch_type w = acos(z); + return { abs(w.imag()), copysign(w.real(), z.imag()) }; } // asin @@ -130,20 +186,8 @@ namespace xsimd XSIMD_INLINE batch, A> asin(const batch, A>& z, requires_arch) noexcept { using batch_type = batch, A>; - using real_batch = typename batch_type::real_batch; - real_batch x = z.real(); - real_batch y = z.imag(); - - batch_type ct(-y, x); - batch_type zz(real_batch(1.) - (x - y) * (x + y), -2 * x * y); - zz = log(ct + sqrt(zz)); - batch_type resg(zz.imag(), -zz.real()); - - return select(y == real_batch(0.), - select(fabs(x) > real_batch(1.), - batch_type(constants::pio2(), real_batch(0.)), - batch_type(asin(x), real_batch(0.))), - resg); + const batch_type w = detail::asin_kernel(z); + return { atan2(z.real(), w.real()), copysign(w.imag(), z.imag()) }; } // asinh @@ -354,15 +398,14 @@ namespace xsimd XSIMD_INLINE batch atan2(batch const& self, batch const& other, requires_arch) noexcept { using batch_type = batch; + const batch_type q = select(self == batch_type(0.), batch_type(0.), abs(self / other)); // atan2(+-0, +-0): 0 / 0 would be NaN #ifdef __FAST_MATH__ - const batch_type q = abs(self / other); const batch_type q_p = abs(other / self); #else - const batch_type q = abs(self / 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); + return select(signnz(other) > batch_type(0.), z, constants::pi() - z) * signnz(self); } // cos @@ -858,21 +901,13 @@ namespace xsimd const batch_type y = detail::tan_eval(xr, test); return y ^ bitofsign(self); } + template XSIMD_INLINE batch, A> tan(batch, A> const& z, requires_arch) noexcept { using batch_type = batch, A>; - using real_batch = typename batch_type::real_batch; - real_batch d = cos(2 * z.real()) + cosh(2 * z.imag()); - real_batch wreal = sin(2 * z.real()) / d; - real_batch wimag = sinh(2 * z.imag()); -#ifdef __FAST_MATH__ - return batch_type(wreal, real_batch(1.)), batch_type(wreal, wimag / d); -#else - batch_type winf(constants::infinity(), constants::infinity()); - batch_type wres = select(isinf(wimag), batch_type(wreal, real_batch(1.)), batch_type(wreal, wimag / d)); - return select(d == real_batch(0.), winf, wres); -#endif + const batch_type w = tanh(batch_type(-z.imag(), z.real())); + return { w.imag(), -w.real() }; } // tanh @@ -978,15 +1013,46 @@ namespace xsimd batch_type r = fma(batch_type(-2.), one / (one + exp(x + x)), one); return select(test, z, r) ^ bts; } + // With z = x + iy, + // tanh(z) = ((1 + tan(y)^2) sinh(x) cosh(x) + i tan(y)) / (1 + (1 + tan(y)^2) sinh(x)^2) + // as in FreeBSD's ctanh: https://github.com/freebsd/freebsd-src/blob/main/lib/msun/src/s_ctanh.c + // It is (sinh(2x) + i sin(2y)) / (cosh(2x) + cos(2y)) with numerator and denominator divided by + // 2 cos(y)^2. For |x| >= 11 (float) or 22 (double) it rounds to + // tanh(z) = sign(x) + i tan(y) / ((1 + tan(y)^2) sinh(x)^2) + // + // In the code, with x = Re z and y = Im z: + // t = tan(y), beta = 1 + t^2, s = sinh(|x|), c = cosh(|x|) = sqrt(1 + s^2) + // tanh(z) = sign(x) beta s c / (1 + beta s^2) + i t / (1 + beta s^2) + // tanh(z) = sign(x) + i t / (beta s^2) for |x| >= 11 (float) or 22 (double) template XSIMD_INLINE batch, A> tanh(const batch, A>& z, requires_arch) noexcept { - using real_batch = typename batch, A>::real_batch; - auto x = z.real(); - auto y = z.imag(); - real_batch two(2); - auto d = cosh(two * x) + cos(two * y); - return { sinh(two * x) / d, sin(two * y) / d }; + using batch_type = batch, A>; + using real_batch = typename batch_type::real_batch; + const real_batch one(1.); + const real_batch x = z.real(); + const real_batch ax = abs(x); + + // t and beta = 1 / cos(y)^2 >= 1. + const real_batch t = tan(z.imag()); + const real_batch beta = fma(t, t, one); + + // s and c. |x| is capped so that s stays finite: with -ffast-math, sinh of a large |x| is + // garbage instead of inf, and 1 / inf is NaN. + const real_batch cap(std::is_same_v ? T(64) : T(512)); + const real_batch s = sinh(select(ax > cap, cap, ax)); + const real_batch c = sqrt(fma(s, s, one)); + + // Main formula. All terms are >= 0 and the denominator is >= 1: no cancellation, no division by 0. + const real_batch inv = one / fma(beta * s, s, one); + const real_batch re = beta * s * c * inv; + const real_batch im = t * inv; + + // Large |x|. There 1 + beta s^2 rounds to beta s^2 and c / s rounds to 1, so re = sign(x) and + // im = t / (beta s^2), computed as (t / beta) (1 / s) (1 / s) because s^2 can overflow. + const auto large = ax >= real_batch(std::is_same_v ? T(11) : T(22)); + const real_batch inv_s = one / s; + return { copysign(select(large, one, re), x), select(large, (t / beta) * inv_s * inv_s, im) }; } } diff --git a/include/xsimd/arch/xsimd_scalar.hpp b/include/xsimd/arch/xsimd_scalar.hpp index 3e2aa13a4..ecfc744f0 100644 --- a/include/xsimd/arch/xsimd_scalar.hpp +++ b/include/xsimd/arch/xsimd_scalar.hpp @@ -701,7 +701,9 @@ namespace xsimd { using T = typename C::value_type; C u = C(1.) + val; - return u == C(1.) ? val : (u.real() <= T(0.) ? log(u) : log(u) * val / (u - C(1.))); + C res = u == C(1.) ? val : (u.real() <= T(0.) ? log(u) : log(u) * val / (u - C(1.))); + // Im(log(1 + z)) has the sign of Im(z), which the operations above lose for an imaginary -0. + return C(res.real(), std::copysign(res.imag(), val.imag())); } } diff --git a/test/test_complex_exponential.cpp b/test/test_complex_exponential.cpp index 8dffdf5c4..9a41a6531 100644 --- a/test/test_complex_exponential.cpp +++ b/test/test_complex_exponential.cpp @@ -213,4 +213,15 @@ TEST_CASE_TEMPLATE("[complex exponential]", B, BATCH_COMPLEX_TYPES) Test.test_sign(); } } + +TEST_CASE_TEMPLATE("[complex exponential issue 1432]", B, BATCH_COMPLEX_TYPES) +{ + using value_type = typename B::value_type; + for (const auto& z : detail::issue_1432_inputs()) + { + const B expected(std::log(value_type(1 + z.real(), z.imag()))); + const B actual = log1p(B(z)); + CHECK_BATCH_EQ(expected, actual); + } +} #endif diff --git a/test/test_complex_trigonometric.cpp b/test/test_complex_trigonometric.cpp index 89fdbd4aa..78959bc74 100644 --- a/test/test_complex_trigonometric.cpp +++ b/test/test_complex_trigonometric.cpp @@ -229,4 +229,40 @@ TEST_CASE_TEMPLATE("[complex trigonometric]", B, BATCH_COMPLEX_TYPES) Test.test_atan(); } } + +// Complex functions from xsimd_common_trigo.hpp on the inputs of issue 1432, compared with std. +TEST_CASE_TEMPLATE("[complex trigonometric issue 1432]", B, BATCH_COMPLEX_TYPES) +{ + const auto inputs = detail::issue_1432_inputs(); + + SUBCASE("asin") + { + CHECK_AGAINST_STD(B, asin, inputs); + } + + SUBCASE("acos") + { + CHECK_AGAINST_STD(B, acos, inputs); + } + + SUBCASE("tan") + { + CHECK_AGAINST_STD(B, tan, inputs); + } + + SUBCASE("asinh") + { + CHECK_AGAINST_STD(B, asinh, inputs); + } + + SUBCASE("acosh") + { + CHECK_AGAINST_STD(B, acosh, inputs); + } + + SUBCASE("tanh") + { + CHECK_AGAINST_STD(B, tanh, inputs); + } +} #endif diff --git a/test/test_trigonometric.cpp b/test/test_trigonometric.cpp index d3e8862a6..4aee1f012 100644 --- a/test/test_trigonometric.cpp +++ b/test/test_trigonometric.cpp @@ -178,6 +178,19 @@ struct trigonometric_test CHECK_BATCH_EQ(ref, out); } } + // atan2 with both arguments zero: 0, +-pi as in std, not NaN + { + const value_type zeros[] = { value_type(0.), value_type(-0.) }; + for (const value_type x : zeros) + { + for (const value_type y : zeros) + { + const batch_type ref(std::atan2(x, y)); + const batch_type out = atan2(batch_type(x), batch_type(y)); + CHECK_BATCH_EQ(ref, out); + } + } + } } }; diff --git a/test/test_utils.hpp b/test/test_utils.hpp index f0d64a7df..d1a4fcaa1 100644 --- a/test/test_utils.hpp +++ b/test/test_utils.hpp @@ -468,6 +468,40 @@ namespace detail } \ } while (0) +// Checks xsimd::f against std::f on every element of inputs, each broadcast to a whole batch of type B. +#define CHECK_AGAINST_STD(B, f, inputs) \ + do \ + { \ + for (const auto& input : inputs) \ + { \ + const B expected(std::f(input)); \ + const B actual = xsimd::f(B(input)); \ + CHECK_BATCH_EQ(expected, actual); \ + } \ + } while (0) + +namespace detail +{ + // Inputs from https://github.com/xtensor-stack/xsimd/issues/1432 + template + std::vector issue_1432_inputs() + { + using value_type = typename B::value_type; + using real_value_type = typename value_type::value_type; + const real_value_type re[] = { -400, -100, -42.4117, -4, -2, -1, -0., 0., 1, 1.5707, 1.57076, 2, 45, 356, 633.708, 700, 10000 }; + const real_value_type im[] = { -500, -356, -100, -45, -10.9958, -3, -1, -0., 0., 1, 2, 5, 800 }; + std::vector inputs; + for (const real_value_type x : re) + { + for (const real_value_type y : im) + { + inputs.emplace_back(x, y); + } + } + return inputs; + } +} + /*********************** * Testing types lists * ***********************/ From ac774c8853bf988fee7658864dc7ce120a9c6a17 Mon Sep 17 00:00:00 2001 From: Fikret Ardal Date: Thu, 12 Nov 2026 12:16:41 -0500 Subject: [PATCH 2/2] Test complex issue 1432 functions on hardcoded cases Check asin, acos, asinh, acosh, tan, tanh and log1p on the inputs of issue 1432 that gave a wrong value or a wrong sign of zero before the fix: signed zeros, inputs where the sign of a zero imaginary part decides the sign of the result, points near the poles of tan and tanh, and large |z|. The expected values are the exact results rounded to double, written in the tests. A result passes if it is within 4 * cond * eps of the expected value, where cond is the condition number |z f'(z) / f(z)|. Infinite values and the sign of zeros must match; with fast-math, infinite values are skipped. Cite the papers whose formulas complex asin and tanh use (Kahan 1987; Hull, Fairgrieve and Tang 1997) instead of the Julia and FreeBSD code. Assisted-By: Claude Opus 5.5 --- .../xsimd/arch/common/xsimd_common_trigo.hpp | 10 +- test/test_complex_exponential.cpp | 26 ++++-- test/test_complex_trigonometric.cpp | 92 +++++++++++++++++-- test/test_utils.hpp | 58 +++++++----- 4 files changed, 142 insertions(+), 44 deletions(-) diff --git a/include/xsimd/arch/common/xsimd_common_trigo.hpp b/include/xsimd/arch/common/xsimd_common_trigo.hpp index 6119c65fb..48d14a6fb 100644 --- a/include/xsimd/arch/common/xsimd_common_trigo.hpp +++ b/include/xsimd/arch/common/xsimd_common_trigo.hpp @@ -51,7 +51,8 @@ namespace xsimd // For z = x + iy, returns Re sqrt(1 - z^2) as the real part and |Im asin(z)| as the imaginary part. // With p = sqrt(1 + z) and q = sqrt(1 - z): // Re sqrt(1 - z^2) = Re(q * p), |Im asin(z)| = asinh(Im(conj(q) p)) - // as in Julia's asin(z::Complex): https://github.com/JuliaLang/julia/blob/master/base/complex.jl + // from W. Kahan, "Branch Cuts for Complex Elementary Functions, or Much Ado About Nothing's Sign Bit", + // in The State of the Art in Numerical Analysis, A. Iserles and M. J. D. Powell (eds.), 1987, pp. 165-211. template XSIMD_INLINE batch, A> asin_kernel(const batch, A>& z) noexcept { @@ -90,8 +91,9 @@ namespace xsimd real_batch v = asinh(sh); // For huge x or y, (1 +- x)^2 + y^2 overflows. There Re sqrt(1 - z^2) = y and - // |Im asin(z)| = log(2 sqrt(x^2 + y^2)) up to rounding, as in FreeBSD's casinh: - // https://github.com/freebsd/freebsd-src/blob/main/lib/msun/src/catrig.c + // |Im asin(z)| = log(2 sqrt(x^2 + y^2)) up to rounding, as in T. E. Hull, T. F. Fairgrieve and + // P. T. P. Tang, "Implementing the Complex Arcsine and Arccosine Functions Using Exception Handling", + // ACM Trans. Math. Software 23(3), 1997, pp. 299-335. const auto y_larger = x < y; const real_batch m = select(y_larger, y, x); // max(|x|, |y|) const auto huge = m >= constants::oneotwoeps(); @@ -1015,7 +1017,7 @@ namespace xsimd } // With z = x + iy, // tanh(z) = ((1 + tan(y)^2) sinh(x) cosh(x) + i tan(y)) / (1 + (1 + tan(y)^2) sinh(x)^2) - // as in FreeBSD's ctanh: https://github.com/freebsd/freebsd-src/blob/main/lib/msun/src/s_ctanh.c + // from W. Kahan, "Branch Cuts for Complex Elementary Functions, or Much Ado About Nothing's Sign Bit" (1987). // It is (sinh(2x) + i sin(2y)) / (cosh(2x) + cos(2y)) with numerator and denominator divided by // 2 cos(y)^2. For |x| >= 11 (float) or 22 (double) it rounds to // tanh(z) = sign(x) + i tan(y) / ((1 + tan(y)^2) sinh(x)^2) diff --git a/test/test_complex_exponential.cpp b/test/test_complex_exponential.cpp index 9a41a6531..43930b64f 100644 --- a/test/test_complex_exponential.cpp +++ b/test/test_complex_exponential.cpp @@ -214,14 +214,26 @@ TEST_CASE_TEMPLATE("[complex exponential]", B, BATCH_COMPLEX_TYPES) } } +// Cases of issue 1432 that the fix changed, as in test_complex_trigonometric.cpp. TEST_CASE_TEMPLATE("[complex exponential issue 1432]", B, BATCH_COMPLEX_TYPES) { - using value_type = typename B::value_type; - for (const auto& z : detail::issue_1432_inputs()) - { - const B expected(std::log(value_type(1 + z.real(), z.imag()))); - const B actual = log1p(B(z)); - CHECK_BATCH_EQ(expected, actual); - } + using xsimd::constants::log_2; + using xsimd::constants::pi; + + std::vector> inputs = { + { 1, -0. }, // imaginary part -0 + { -2, -0. }, // the -0 sets the sign of Im log1p(z) + }; + std::vector> expected = { + { log_2(), -0. }, + { 0., -pi() }, + }; + std::vector cond = { 1, 1 }; +#ifndef __FAST_MATH__ // no infinities with fast-math + inputs.push_back({ -1, 0. }); + expected.push_back({ -std::numeric_limits::infinity(), 0. }); + cond.push_back(1); +#endif + CHECK_COMPLEX_CASES(B, log1p, inputs, expected, cond); } #endif diff --git a/test/test_complex_trigonometric.cpp b/test/test_complex_trigonometric.cpp index 78959bc74..152d561d0 100644 --- a/test/test_complex_trigonometric.cpp +++ b/test/test_complex_trigonometric.cpp @@ -230,39 +230,115 @@ TEST_CASE_TEMPLATE("[complex trigonometric]", B, BATCH_COMPLEX_TYPES) } } -// Complex functions from xsimd_common_trigo.hpp on the inputs of issue 1432, compared with std. +// Cases of issue 1432 that the fix changed: each gave a wrong value or a wrong sign of zero before it. The expected +// values are the exact f(z) rounded to double; cond is max(1, |z f'(z) / f(z)|) (see detail::check_complex_cases). TEST_CASE_TEMPLATE("[complex trigonometric issue 1432]", B, BATCH_COMPLEX_TYPES) { - const auto inputs = detail::issue_1432_inputs(); + using xsimd::constants::pi; + using xsimd::constants::pio2; SUBCASE("asin") { - CHECK_AGAINST_STD(B, asin, inputs); + const std::vector> inputs = { + { -2, -0. }, // the -0 sets the sign of Im asin(z) + { -0., -0. }, // signed zeros + { 1, 800 }, // large imaginary part + { 10000, 1 }, // large real part + }; + const std::vector> expected = { + { -pio2(), -1.3169578969248168 }, + { -0., -0. }, + { 0.0012499983724006398, 7.377760080100202 }, + { 1.57069632679473, 9.903487555036127 }, + }; + const std::vector cond = { 1, 1, 1, 1 }; + CHECK_COMPLEX_CASES(B, asin, inputs, expected, cond); } SUBCASE("acos") { - CHECK_AGAINST_STD(B, acos, inputs); + const std::vector> inputs = { + { -2, -0. }, // the -0 sets the sign of Im acos(z) + { 0., -0. }, // signed zeros + { 1, 800 }, // large imaginary part + { 10000, 1 }, // large real part + }; + const std::vector> expected = { + { pi(), 1.3169578969248168 }, + { pio2(), 0. }, + { 1.569546328422496, -7.377760080100202 }, + { 0.00010000000016666667, -9.903487555036127 }, + }; + const std::vector cond = { 1, 1, 1, 1 }; + CHECK_COMPLEX_CASES(B, acos, inputs, expected, cond); } SUBCASE("tan") { - CHECK_AGAINST_STD(B, tan, inputs); + const std::vector> inputs = { + { 1.5707f, 0. }, // near the pole pi/2 = 1.570796 + { 2, -500 }, // large imaginary part + }; + const std::vector> expected = { + { 10386.654119664716, 0. }, + { -0., -1. }, + }; + const std::vector cond = { 1.63e4, 1 }; + CHECK_COMPLEX_CASES(B, tan, inputs, expected, cond); } SUBCASE("asinh") { - CHECK_AGAINST_STD(B, asinh, inputs); + const std::vector> inputs = { + { -0., -3 }, // the -0 sets the sign of Re asinh(z) + { -0., -0. }, // signed zeros + { 2, 0. }, // real axis + { 1, -500 }, // large imaginary part + { 10000, 1 }, // large real part + }; + const std::vector> expected = { + { -1.762747174039086, -pio2() }, + { -0., -0. }, + { 1.4436354751788103, 0. }, + { 6.907756278988637, -1.5687963254615769 }, + { 9.903487560036128, 9.999999916666668e-05 }, + }; + const std::vector cond = { 1, 1, 1, 1, 1 }; + CHECK_COMPLEX_CASES(B, asinh, inputs, expected, cond); } SUBCASE("acosh") { - CHECK_AGAINST_STD(B, acosh, inputs); + const std::vector> inputs = { + { -4, -3 }, // x < 0, y < 0 + { -2, 0. }, // the +0 sets the sign of Im acosh(z) + { 1, -0. }, // signed zeros + { 1, 800 }, // large imaginary part + { 10000, 1 }, // large real part + }; + const std::vector> expected = { + { 2.2999140408792695, -2.488413180146375 }, + { 1.3169578969248168, pi() }, + { 0., -0. }, + { 7.377760080100202, 1.569546328422496 }, + { 9.903487555036127, 0.00010000000016666667 }, + }; + const std::vector cond = { 1, 1, 1, 1, 1 }; + CHECK_COMPLEX_CASES(B, acosh, inputs, expected, cond); } SUBCASE("tanh") { - CHECK_AGAINST_STD(B, tanh, inputs); + const std::vector> inputs = { + { -400, 2 }, // large real part + { 0., -10.9958f }, // near the pole -7 pi/2 i = -10.995574 i + }; + const std::vector> expected = { + { -1., -0. }, + { 0., 4430.056602996832 }, + }; + const std::vector cond = { 1, 4.87e4 }; + CHECK_COMPLEX_CASES(B, tanh, inputs, expected, cond); } } #endif diff --git a/test/test_utils.hpp b/test/test_utils.hpp index d1a4fcaa1..2a5d08eed 100644 --- a/test/test_utils.hpp +++ b/test/test_utils.hpp @@ -468,40 +468,48 @@ namespace detail } \ } while (0) -// Checks xsimd::f against std::f on every element of inputs, each broadcast to a whole batch of type B. -#define CHECK_AGAINST_STD(B, f, inputs) \ - do \ - { \ - for (const auto& input : inputs) \ - { \ - const B expected(std::f(input)); \ - const B actual = xsimd::f(B(input)); \ - CHECK_BATCH_EQ(expected, actual); \ - } \ - } while (0) - namespace detail { - // Inputs from https://github.com/xtensor-stack/xsimd/issues/1432 - template - std::vector issue_1432_inputs() - { - using value_type = typename B::value_type; - using real_value_type = typename value_type::value_type; - const real_value_type re[] = { -400, -100, -42.4117, -4, -2, -1, -0., 0., 1, 1.5707, 1.57076, 2, 45, 356, 633.708, 700, 10000 }; - const real_value_type im[] = { -500, -356, -100, -45, -10.9958, -3, -1, -0., 0., 1, 2, 5, 800 }; - std::vector inputs; - for (const real_value_type x : re) + template + void check_complex_cases(const char* name, const std::vector>& inputs, + const std::vector>& expected, const std::vector& cond, F f) + { + using T = typename B::value_type::value_type; + const bool sizes_match = expected.size() == inputs.size() && cond.size() == inputs.size(); + if (!sizes_match) + { + FAIL_CHECK(name << ": inputs, expected and cond must have the same size"); + return; + } + for (size_t k = 0; k < inputs.size(); ++k) { - for (const real_value_type y : im) + const std::complex z = inputs[k]; + const std::complex e = expected[k]; + const auto r = f(B(typename B::value_type(T(z.real()), T(z.imag())))).get(0); + const std::complex out(r.real(), r.imag()); + const double allowed = 4 * cond[k] * std::numeric_limits::epsilon() * std::abs(e); +#ifdef __FAST_MATH__ + const bool ok = std::abs(out - e) <= allowed; // no infinite values with fast-math +#else + bool ok = std::isinf(e.real()) || std::isinf(e.imag()) ? out == e : std::abs(out - e) <= allowed; + // zeros must have the right sign (fast-math ignores it) + ok = ok && (e.real() != 0 || std::signbit(out.real()) == std::signbit(e.real())); + ok = ok && (e.imag() != 0 || std::signbit(out.imag()) == std::signbit(e.imag())); +#endif + if (!ok) { - inputs.emplace_back(x, y); + std::ostringstream msg; + msg << std::setprecision(std::numeric_limits::max_digits10) << name << z << " = " << out << ", expected " << e; + FAIL_CHECK(msg.str()); } } - return inputs; } } +// Checks xsimd::f at inputs against expected, within cond (see detail::check_complex_cases), for the value type of B. +#define CHECK_COMPLEX_CASES(B, f, inputs, expected, cond) \ + ::detail::check_complex_cases(#f, inputs, expected, cond, [](const B& z) { return xsimd::f(z); }) + /*********************** * Testing types lists * ***********************/