diff --git a/include/xtensor/core/xmath.hpp b/include/xtensor/core/xmath.hpp index e40536651..2ef3074b7 100644 --- a/include/xtensor/core/xmath.hpp +++ b/include/xtensor/core/xmath.hpp @@ -339,6 +339,22 @@ namespace xt XTENSOR_UNARY_MATH_FUNCTOR(isinf); XTENSOR_UNARY_MATH_FUNCTOR(isnan); XTENSOR_UNARY_MATH_FUNCTOR(conj); + + struct sinc_fun + { + template + auto operator()(const T& x) const + { + using math::sin; + using result_type = decltype(sin(x)); + if (x == result_type(0)) + { + return result_type(1); + } + const auto angle = numeric_constants::PI * result_type(x); + return sin(angle) / angle; + } + }; } #undef XTENSOR_UNARY_MATH_FUNCTOR @@ -1259,6 +1275,21 @@ namespace xt return detail::make_xfunction(std::forward(e)); } + /** + * @ingroup trigo_functions + * @brief Normalized sinc function. + * + * Computes \f$\sin(\pi x)/(\pi x)\f$ element-wise, with value 1 at zero. + * This is the normalized convention used by NumPy. + * @param e an \ref xexpression + * @return an \ref xfunction + */ + template + inline auto sinc(E&& e) noexcept -> detail::xfunction_type_t + { + return detail::make_xfunction(std::forward(e)); + } + /** * @ingroup trigo_functions * @brief Cosine function. diff --git a/test/test_xmath.cpp b/test/test_xmath.cpp index 982e71a1a..0f8174754 100644 --- a/test/test_xmath.cpp +++ b/test/test_xmath.cpp @@ -12,9 +12,11 @@ #include "xtensor/containers/xadapt.hpp" #include "xtensor/containers/xarray.hpp" +#include "xtensor/containers/xscalar.hpp" #include "xtensor/core/xmath.hpp" #include "xtensor/generators/xrandom.hpp" #include "xtensor/optional/xoptional_assembly.hpp" +#include "xtensor/views/xview.hpp" #include "test_common_macros.hpp" @@ -540,6 +542,109 @@ namespace xt EXPECT_EQ(sin(a)(0, 0), std::sin(a(0, 0))); } + TEST(xmath, sinc) + { + xarray input = {{-1.5, -0.5, -1e-8}, {0.0, 0.5, 1.5}}; + const auto original = input; + const xarray expected = { + {-0.2122065907891938, 0.6366197723675814, 0.9999999999999999}, + {1.0, 0.6366197723675814, -0.2122065907891938} + }; + const xarray result = sinc(input); + EXPECT_EQ(result.shape(), input.shape()); + EXPECT_TRUE(allclose(result, expected, 1e-14, 1e-15)); + EXPECT_EQ(input, original); + + // Rounded NumPy reference values reported in issue1608. + xarray reference_input = {-4.0, -3.8, -3.6, -0.2, 0.0, 0.2, 3.6, 3.8, 4.0}; + xarray reference = {-3.89804309e-17, -4.92362781e-02, -8.40918587e-02, + 9.35489284e-01, 1.0, 9.35489284e-01, + -8.40918587e-02, -4.92362781e-02, -3.89804309e-17}; + EXPECT_TRUE(allclose(sinc(reference_input), reference, 1e-8, 1e-9)); + + xarray values = {-0.0f, 0.5f}; + auto expression = sinc(xt::view(values, xt::all())); + static_assert(std::is_same_v); + EXPECT_EQ(expression(0), 1.0f); + values(0) = -0.5f; + EXPECT_LE(std::abs(expression(0) - 0.63661977f), 1e-7f); + + xarray integers = {0, 1}; + static_assert(std::is_same_v); + EXPECT_EQ(sinc(integers)(0), 1.0); + EXPECT_LE(std::abs(sinc(integers)(1)), 1e-15); + } + + TEST(xmath, sinc_expressions) + { + xscalar scalar_integer(0); + xscalar scalar_double(0.5); + xscalar scalar_float(0.0f); + EXPECT_EQ(sinc(scalar_integer)(), 1.0); + EXPECT_LE(std::abs(sinc(scalar_double)() - 2.0 / numeric_constants::PI), 1e-15); + static_assert(std::is_same_v); + static_assert(std::is_same_v); + + xarray input = {{-1.5, -0.5}, {0.5, 1.5}}; + const auto original = input; + auto expression = sinc(input + 0.0); + EXPECT_EQ(expression.shape(), input.shape()); + EXPECT_TRUE(allclose(expression, sinc(-input), 1e-14, 1e-15)); + input(0, 0) = 0.0; + EXPECT_EQ(expression(0, 0), 1.0); + EXPECT_LE(std::abs(expression(0, 1) - 2.0 / numeric_constants::PI), 1e-15); + + xarray small = {-1e-8f, -0.0f, 0.0f, 1e-8f}; + const xarray near_zero = sinc(small); + EXPECT_TRUE(allclose(near_zero, ones({4}), 1e-6f, 1e-7f)); + EXPECT_EQ(small(1), -0.0f); + EXPECT_EQ(original(0, 0), -1.5); + + xarray empty = xarray::from_shape({0, 2}); + const xarray empty_result = sinc(empty); + EXPECT_EQ(empty_result.shape(), empty.shape()); + EXPECT_EQ(empty_result.size(), 0u); + } + + TEST(xmath, sinc_complex) + { + using complex_type = std::complex; + xarray input = {complex_type(0.0, 0.0), complex_type(0.0, 0.5), + complex_type(0.5, 0.5)}; + const auto original = input; + const xarray result = sinc(input); + static_assert(std::is_same_v); + EXPECT_EQ(result(0), complex_type(1.0, 0.0)); + EXPECT_LE(std::abs(result(1).real() - std::sinh(numeric_constants::PI_2) + / numeric_constants::PI_2), 1e-14); + EXPECT_LE(std::abs(result(1).imag()), 1e-15); + const double expected = std::cosh(numeric_constants::PI_2) + / numeric_constants::PI; + EXPECT_LE(std::abs(result(2).real() - expected), 1e-14); + EXPECT_LE(std::abs(result(2).imag() + expected), 1e-14); + EXPECT_EQ(input, original); + } + + TEST(xmath, sinc_optional) + { + using opt_type = xoptional_assembly, xarray>; + opt_type input = {0.0, 0.5, xtl::missing()}; + auto expression = sinc(input); + // Direct access exercises reference-backed optional proxies as well as + // the split value/missing-mask path used by assignment. + EXPECT_TRUE(expression(0).has_value()); + EXPECT_EQ(expression(0).value(), 1.0); + EXPECT_TRUE(expression(1).has_value()); + EXPECT_LE(std::abs(expression(1).value() - 2.0 / numeric_constants::PI), 1e-15); + EXPECT_FALSE(expression(2).has_value()); + const opt_type result = expression; + EXPECT_EQ(result(0).value(), 1.0); + EXPECT_LE(std::abs(result(1).value() - 2.0 / numeric_constants::PI), 1e-15); + EXPECT_FALSE(result(2).has_value()); + input(0) = 0.5; + EXPECT_LE(std::abs(expression(0).value() - 2.0 / numeric_constants::PI), 1e-15); + } + TEST(xmath, cos) { shape_type shape = {3, 2};