diff --git a/include/xtensor/core/xfunction.hpp b/include/xtensor/core/xfunction.hpp index 99144b377..0c94b4fe2 100644 --- a/include/xtensor/core/xfunction.hpp +++ b/include/xtensor/core/xfunction.hpp @@ -612,7 +612,7 @@ namespace xt { XTENSOR_TRY(check_index(shape(), args...)); XTENSOR_CHECK_DIMENSION(shape(), args...); - return m_f(e(args...)...); + return m_f(e(static_cast(args)...)...); }, m_e ); diff --git a/include/xtensor/core/xmath.hpp b/include/xtensor/core/xmath.hpp index 6b86c4b74..07a9b09a0 100644 --- a/include/xtensor/core/xmath.hpp +++ b/include/xtensor/core/xmath.hpp @@ -18,6 +18,7 @@ #include #include #include +#include #include #include @@ -38,18 +39,18 @@ namespace xt template struct numeric_constants { - static constexpr T PI = 3.141592653589793238463; - static constexpr T PI_2 = 1.57079632679489661923; - static constexpr T PI_4 = 0.785398163397448309616; - static constexpr T D_1_PI = 0.318309886183790671538; - static constexpr T D_2_PI = 0.636619772367581343076; - static constexpr T D_2_SQRTPI = 1.12837916709551257390; - static constexpr T SQRT2 = 1.41421356237309504880; - static constexpr T SQRT1_2 = 0.707106781186547524401; - static constexpr T E = 2.71828182845904523536; - static constexpr T LOG2E = 1.44269504088896340736; - static constexpr T LOG10E = 0.434294481903251827651; - static constexpr T LN2 = 0.693147180559945309417; + static constexpr T PI = std::numbers::pi_v; + static constexpr T PI_2 = std::numbers::pi_v / 2; + static constexpr T PI_4 = std::numbers::pi_v / 4; + static constexpr T D_1_PI = std::numbers::inv_pi_v; + static constexpr T D_2_PI = 2 * std::numbers::inv_pi_v; + static constexpr T D_2_SQRTPI = 2 * std::numbers::inv_sqrtpi_v; + static constexpr T SQRT2 = std::numbers::sqrt2_v; + static constexpr T SQRT1_2 = std::numbers::sqrt2_v / 2; + static constexpr T E = std::numbers::e_v; + static constexpr T LOG2E = std::numbers::log2e_v; + static constexpr T LOG10E = std::numbers::log10e_v; + static constexpr T LN2 = std::numbers::ln2_v; }; /*********** diff --git a/include/xtensor/io/xnpy.hpp b/include/xtensor/io/xnpy.hpp index ef517f068..b902f894d 100644 --- a/include/xtensor/io/xnpy.hpp +++ b/include/xtensor/io/xnpy.hpp @@ -479,9 +479,11 @@ namespace xt char header_len_le16[2]; istream.read(header_len_le16, 2); - uint16_t header_b0 = static_cast(static_cast(header_len_le16[0])); - uint16_t header_b1 = static_cast(static_cast(header_len_le16[1])) << 8; - uint16_t header_length = header_b0 | header_b1; + // The length is stored on 16 bits. It is assembled in a uint32_t like in + // read_header_2_0, because arithmetic on uint16_t is done in int. + uint32_t header_b0 = static_cast(static_cast(header_len_le16[0])); + uint32_t header_b1 = static_cast(static_cast(header_len_le16[1])) << 8; + uint32_t header_length = header_b0 | header_b1; if ((magic_string_length + 2 + 2 + header_length) % 16 != 0) { diff --git a/include/xtensor/misc/xfft.hpp b/include/xtensor/misc/xfft.hpp index caf37145e..8d40fa6ba 100644 --- a/include/xtensor/misc/xfft.hpp +++ b/include/xtensor/misc/xfft.hpp @@ -61,8 +61,10 @@ namespace xt auto odd = radix2(xt::view(ev, xt::range(1, _, 2))); #endif - auto range = xt::arange(static_cast(N) / 2); - auto exp = xt::exp(static_cast(-2i) * pi * range / N); + // N is a power of 2, so it is exactly representable as a double. + const auto N_dbl = static_cast(N); + auto range = xt::arange(N_dbl / 2); + auto exp = xt::exp(static_cast(-2i) * pi * range / N_dbl); auto t = exp * odd; auto first_half = even + t; auto second_half = even - t; @@ -90,7 +92,10 @@ namespace xt xt::xtensor i = xt::pow(xt::linspace(0, n - 1, n), 2); i %= (n * 2); - auto angles = xt::eval(precision{3.141592653589793238463} * i / n); + // The angles are computed in precision, the integer operands are converted to it. + auto angles = xt::eval( + xt::numeric_constants::PI * xt::cast(i) / static_cast(n) + ); auto j = std::complex(0, 1); exp_table = xt::exp(-angles * j); @@ -112,7 +117,8 @@ namespace xt auto spectrum_k = xv * yv; auto complex_args = xt::conj(spectrum_k); auto fft_res = radix2(complex_args); - auto cv = xt::conj(fft_res) / m; + // m is a power of 2, so it is exactly representable in floating point. + auto cv = xt::conj(fft_res) / static_cast(m); return xt::eval(xt::view(cv, xt::range(0, n)) * exp_table); } @@ -226,7 +232,8 @@ namespace xt auto outvec = ifft(xv, axis); // Scaling (because this FFT implementation omits it) - outvec = outvec / n; + using precision = typename std::decay_t::value_type::value_type; + outvec = outvec / static_cast(n); return outvec; } diff --git a/include/xtensor/misc/xmanipulation.hpp b/include/xtensor/misc/xmanipulation.hpp index db4125e01..6ac6f1067 100644 --- a/include/xtensor/misc/xmanipulation.hpp +++ b/include/xtensor/misc/xmanipulation.hpp @@ -319,8 +319,8 @@ namespace xt check_axis_in_dim(axis1, dim, "Parameter axis1"); check_axis_in_dim(axis2, dim, "Parameter axis2"); - using strides_t = get_strides_t::shape_type>; - return transpose(std::forward(e), detail::swapaxes_perm(dim, axis1, axis2)); + using permutation_t = xindex_type_t::shape_type>; + return transpose(std::forward(e), detail::swapaxes_perm(dim, axis1, axis2)); } /***************************** @@ -373,8 +373,8 @@ namespace xt check_axis_in_dim(src, dim, "Parameter src"); check_axis_in_dim(dest, dim, "Parameter dest"); - using strides_t = get_strides_t::shape_type>; - return xt::transpose(std::forward(e), detail::moveaxis_perm(e.dimension(), src, dest)); + using permutation_t = xindex_type_t::shape_type>; + return xt::transpose(std::forward(e), detail::moveaxis_perm(e.dimension(), src, dest)); } /************************************ diff --git a/include/xtensor/misc/xsort.hpp b/include/xtensor/misc/xsort.hpp index 196ec45cd..d90026530 100644 --- a/include/xtensor/misc/xsort.hpp +++ b/include/xtensor/misc/xsort.hpp @@ -77,7 +77,10 @@ namespace xt { XTENSOR_ASSERT(ev.dimension() >= 2); - const std::size_t n_iters = leading_axis_n_iters(ev); + // Pointer arithmetic needs a signed count. n_iters * secondary_stride is the + // size of the container (or the stride is 0 and n_iters does not matter), + // so n_iters fits in a ptrdiff_t. + const auto n_iters = static_cast(leading_axis_n_iters(ev)); const std::ptrdiff_t secondary_stride = get_secondary_stride(ev); const auto begin = ev.data(); @@ -94,7 +97,8 @@ namespace xt XTENSOR_ASSERT(e1.dimension() >= 2); XTENSOR_ASSERT(e1.dimension() == e2.dimension()); - const std::size_t n_iters = leading_axis_n_iters(e1); + // See above, n_iters fits in a ptrdiff_t. + const auto n_iters = static_cast(leading_axis_n_iters(e1)); const std::ptrdiff_t secondary_stride1 = get_secondary_stride(e1); const std::ptrdiff_t secondary_stride2 = get_secondary_stride(e2); XTENSOR_ASSERT(secondary_stride1 == secondary_stride2); @@ -306,29 +310,25 @@ namespace xt (void) idx_end; // TODO(C++17) [[maybe_unused]] only used in assertion. std::iota(idx_begin, idx_end, 0); + // The indices are in [0, std::distance(data_begin, data_end)), so they fit in + // the difference type of the data iterator. + using difference_type = typename std::iterator_traits::difference_type; + auto index_comp = [&](const auto i, const auto j) + { + return comp( + *(data_begin + static_cast(i)), + *(data_begin + static_cast(j)) + ); + }; switch (method) { case (sorting_method::quick): { - std::sort( - idx_begin, - idx_end, - [&](const auto i, const auto j) - { - return comp(*(data_begin + i), *(data_begin + j)); - } - ); + std::sort(idx_begin, idx_end, index_comp); } case (sorting_method::stable): { - std::stable_sort( - idx_begin, - idx_end, - [&](const auto i, const auto j) - { - return comp(*(data_begin + i), *(data_begin + j)); - } - ); + std::stable_sort(idx_begin, idx_end, index_comp); } } } @@ -820,7 +820,7 @@ namespace xt template inline auto select_indices(const S& shape, const I& indices, std::size_t axis) { - using index_type = get_strides_t; + using index_type = xindex_type_t; auto out = std::vector(); select_indices_impl(shape, indices, axis, 0, xtl::make_sequence(shape.size()), out); return out; @@ -832,7 +832,7 @@ namespace xt inline auto fancy_indexing(E&& e, const I& indices, std::ptrdiff_t axis) { const std::size_t ax = normalize_axis(e.dimension(), axis); - using shape_t = get_strides_t::shape_type>; + using shape_t = xindex_type_t::shape_type>; auto shape = xtl::forward_sequence(e.shape()); shape[ax] = indices.size(); return reshape_view( @@ -907,7 +907,7 @@ namespace xt XTENSOR_ASSERT(0. <= beta); XTENSOR_ASSERT(beta <= 1.); - using tmp_shape_t = get_strides_t::shape_type>; + using tmp_shape_t = xindex_type_t::shape_type>; using id_t = typename tmp_shape_t::value_type; const std::size_t ax = normalize_axis(e.dimension(), axis); @@ -915,8 +915,8 @@ namespace xt auto kth_gamma = detail::quantile_kth_gamma(n, probas, alpha, beta); // Select relevant values for computing interpolating quantiles - auto e_partition = xt::partition(std::forward(e), kth_gamma.first, ax); - auto e_kth = detail::fancy_indexing(std::move(e_partition), std::move(kth_gamma.first), ax); + auto e_partition = xt::partition(std::forward(e), kth_gamma.first, axis); + auto e_kth = detail::fancy_indexing(std::move(e_partition), std::move(kth_gamma.first), axis); // Reshape interpolation coefficients auto gm1_g_shape = xtl::make_sequence(e.dimension(), 1); @@ -931,8 +931,9 @@ namespace xt e_kth_g_shape[ax] = 2; e_kth_g_shape[ax + 1] /= 2; auto quantiles = xt::sum(reshape_view(std::move(e_kth_g), std::move(e_kth_g_shape)), ax); - // Cannot do a transpose on a non-strided expression so we have to eval - return moveaxis(eval(std::move(quantiles)), ax, 0); + // Cannot do a transpose on a non-strided expression so we have to eval. + // quantiles has as many dimensions as e, so axis designates the same axis. + return moveaxis(eval(std::move(quantiles)), axis, 0); } // Static proba array overload diff --git a/include/xtensor/views/xindex_view.hpp b/include/xtensor/views/xindex_view.hpp index 6701d3bf3..a07fa1664 100644 --- a/include/xtensor/views/xindex_view.hpp +++ b/include/xtensor/views/xindex_view.hpp @@ -511,7 +511,9 @@ namespace xt template inline auto xindex_view::element(It first, It /*last*/) -> reference { - return m_e[m_indices[(*first)]]; + // The index may be signed (strided views unravel flat indices into signed + // indices) but it is never negative. + return m_e[m_indices[static_cast(*first)]]; } /** @@ -524,7 +526,9 @@ namespace xt template inline auto xindex_view::element(It first, It /*last*/) const -> const_reference { - return m_e[m_indices[(*first)]]; + // The index may be signed (strided views unravel flat indices into signed + // indices) but it is never negative. + return m_e[m_indices[static_cast(*first)]]; } /** diff --git a/test/test_xfft.cpp b/test/test_xfft.cpp index a1b050ce4..6608a9baf 100644 --- a/test/test_xfft.cpp +++ b/test/test_xfft.cpp @@ -9,10 +9,12 @@ namespace xt { size_t k = 2; size_t n = 8192; - size_t A = 10; + float A = 10; auto x = xt::linspace(0, static_cast(n - 1), n); - xt::xarray y = A * xt::sin(2 * xt::numeric_constants::PI * x * k / n); - auto res = xt::fft::fft(y) / (n / 2); + const auto kf = static_cast(k); + const auto nf = static_cast(n); + xt::xarray y = A * xt::sin(2 * xt::numeric_constants::PI * x * kf / nf); + auto res = xt::fft::fft(y) / static_cast(n / 2); REQUIRE(A == doctest::Approx(std::abs(res(k))).epsilon(.0001)); } @@ -20,10 +22,12 @@ namespace xt { size_t k = 2; size_t n = 8; - size_t A = 10; + float A = 10; auto x = xt::linspace(0, static_cast(n - 1), n); - xt::xarray y = A * xt::sin(2 * xt::numeric_constants::PI * x * k / n); - auto res = xt::fft::ifft(y) / (n / 2); + const auto kf = static_cast(k); + const auto nf = static_cast(n); + xt::xarray y = A * xt::sin(2 * xt::numeric_constants::PI * x * kf / nf); + auto res = xt::fft::ifft(y) / static_cast(n / 2); REQUIRE(A == doctest::Approx(std::abs(res(k))).epsilon(.0001)); } @@ -54,12 +58,14 @@ namespace xt { size_t k = 2; size_t n = 10; - size_t A = 1; + float A = 1; size_t dim = 10; - auto x = xt::linspace(0, n - 1, n) * xt::ones({dim, n}); - xt::xarray y = A * xt::sin(2 * xt::numeric_constants::PI * x * k / n); + auto x = xt::linspace(0, static_cast(n - 1), n) * xt::ones({dim, n}); + const auto kf = static_cast(k); + const auto nf = static_cast(n); + xt::xarray y = A * xt::sin(2 * xt::numeric_constants::PI * x * kf / nf); y = xt::transpose(y); - auto res = xt::fft::fft(y, 0) / (n / 2.0); + auto res = xt::fft::fft(y, 0) / (static_cast(n) / 2); REQUIRE(A == doctest::Approx(std::abs(res(k, 0))).epsilon(.0001)); REQUIRE(A == doctest::Approx(std::abs(res(k, 1))).epsilon(.0001)); } @@ -68,11 +74,13 @@ namespace xt { size_t k = 2; size_t n = 15; - size_t A = 1; + float A = 1; size_t dim = 2; - auto x = xt::linspace(0, n - 1, n) * xt::ones({dim, n}); - xt::xarray y = A * xt::sin(2 * xt::numeric_constants::PI * x * k / n); - auto res = xt::fft::fft(y) / (n / 2.0); + auto x = xt::linspace(0, static_cast(n - 1), n) * xt::ones({dim, n}); + const auto kf = static_cast(k); + const auto nf = static_cast(n); + xt::xarray y = A * xt::sin(2 * xt::numeric_constants::PI * x * kf / nf); + auto res = xt::fft::fft(y) / (static_cast(n) / 2); REQUIRE(A == doctest::Approx(std::abs(res(0, k))).epsilon(.0001)); REQUIRE(A == doctest::Approx(std::abs(res(1, k))).epsilon(.0001)); }