Fix conversion warnings in the library and tests

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<size_type>(args)...)...);
             },
             m_e
         );

include/xtensor/core/xmath.hpp

@@ -18,6 +18,7 @@
 #include <array>
 #include <cmath>
 #include <complex>
+#include <numbers>
 #include <type_traits>
 
 #include <xtl/xcomplex.hpp>
@@ -38,18 +39,18 @@ namespace xt
     template <class T = double>
     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<T>;
+        static constexpr T PI_2 = std::numbers::pi_v<T> / 2;
+        static constexpr T PI_4 = std::numbers::pi_v<T> / 4;
+        static constexpr T D_1_PI = std::numbers::inv_pi_v<T>;
+        static constexpr T D_2_PI = 2 * std::numbers::inv_pi_v<T>;
+        static constexpr T D_2_SQRTPI = 2 * std::numbers::inv_sqrtpi_v<T>;
+        static constexpr T SQRT2 = std::numbers::sqrt2_v<T>;
+        static constexpr T SQRT1_2 = std::numbers::sqrt2_v<T> / 2;
+        static constexpr T E = std::numbers::e_v<T>;
+        static constexpr T LOG2E = std::numbers::log2e_v<T>;
+        static constexpr T LOG10E = std::numbers::log10e_v<T>;
+        static constexpr T LN2 = std::numbers::ln2_v<T>;
     };
 
     /***********

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<uint16_t>(static_cast<uint8_t>(header_len_le16[0]));
-            uint16_t header_b1 = static_cast<uint16_t>(static_cast<uint8_t>(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<uint32_t>(static_cast<uint8_t>(header_len_le16[0]));
+            uint32_t header_b1 = static_cast<uint32_t>(static_cast<uint8_t>(header_len_le16[1])) << 8;
+            uint32_t header_length = header_b0 | header_b1;
 
             if ((magic_string_length + 2 + 2 + header_length) % 16 != 0)
             {

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<double>(static_cast<double>(N) / 2);
-                    auto exp = xt::exp(static_cast<value_type>(-2i) * pi * range / N);
+                    // N is a power of 2, so it is exactly representable as a double.
+                    const auto N_dbl = static_cast<double>(N);
+                    auto range = xt::arange<double>(N_dbl / 2);
+                    auto exp = xt::exp(static_cast<value_type>(-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<std::size_t, 1> i = xt::pow(xt::linspace<std::size_t>(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<precision>::PI * xt::cast<precision>(i) / static_cast<precision>(n)
+                );
                 auto j = std::complex<precision>(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<precision>(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<decltype(outvec)>::value_type::value_type;
+            outvec = outvec / static_cast<precision>(n);
 
             return outvec;
         }

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<typename std::decay_t<E>::shape_type>;
-        return transpose(std::forward<E>(e), detail::swapaxes_perm<strides_t>(dim, axis1, axis2));
+        using permutation_t = xindex_type_t<typename std::decay_t<E>::shape_type>;
+        return transpose(std::forward<E>(e), detail::swapaxes_perm<permutation_t>(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<typename std::decay_t<E>::shape_type>;
-        return xt::transpose(std::forward<E>(e), detail::moveaxis_perm<strides_t>(e.dimension(), src, dest));
+        using permutation_t = xindex_type_t<typename std::decay_t<E>::shape_type>;
+        return xt::transpose(std::forward<E>(e), detail::moveaxis_perm<permutation_t>(e.dimension(), src, dest));
     }
 
     /************************************

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<std::ptrdiff_t>(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<std::ptrdiff_t>(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<ConstRandomIt>::difference_type;
+            auto index_comp = [&](const auto i, const auto j)
+            {
+                return comp(
+                    *(data_begin + static_cast<difference_type>(i)),
+                    *(data_begin + static_cast<difference_type>(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 <class S, class I>
         inline auto select_indices(const S& shape, const I& indices, std::size_t axis)
         {
-            using index_type = get_strides_t<S>;
+            using index_type = xindex_type_t<S>;
             auto out = std::vector<index_type>();
             select_indices_impl(shape, indices, axis, 0, xtl::make_sequence<index_type>(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<typename std::decay_t<E>::shape_type>;
+            using shape_t = xindex_type_t<typename std::decay_t<E>::shape_type>;
             auto shape = xtl::forward_sequence<shape_t, decltype(e.shape())>(e.shape());
             shape[ax] = indices.size();
             return reshape_view(
@@ -907,16 +907,16 @@ namespace xt
         XTENSOR_ASSERT(0. <= beta);
         XTENSOR_ASSERT(beta <= 1.);
 
-        using tmp_shape_t = get_strides_t<typename std::decay_t<E>::shape_type>;
+        using tmp_shape_t = xindex_type_t<typename std::decay_t<E>::shape_type>;
         using id_t = typename tmp_shape_t::value_type;
 
         const std::size_t ax = normalize_axis(e.dimension(), axis);
         const std::size_t n = e.shape()[ax];
         auto kth_gamma = detail::quantile_kth_gamma<T, id_t, P>(n, probas, alpha, beta);
 
         // Select relevant values for computing interpolating quantiles
-        auto e_partition = xt::partition(std::forward<E>(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>(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<tmp_shape_t>(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

include/xtensor/views/xindex_view.hpp

@@ -511,7 +511,9 @@ namespace xt
     template <class It>
     inline auto xindex_view<CT, I>::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<size_type>(*first)]];
     }
 
     /**
@@ -524,7 +526,9 @@ namespace xt
     template <class It>
     inline auto xindex_view<CT, I>::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<size_type>(*first)]];
     }
 
     /**

test/test_xfft.cpp

@@ -9,21 +9,25 @@ namespace xt
     {
         size_t k = 2;
         size_t n = 8192;
-        size_t A = 10;
+        float A = 10;
         auto x = xt::linspace<float>(0, static_cast<float>(n - 1), n);
-        xt::xarray<float> y = A * xt::sin(2 * xt::numeric_constants<float>::PI * x * k / n);
-        auto res = xt::fft::fft(y) / (n / 2);
+        const auto kf = static_cast<float>(k);
+        const auto nf = static_cast<float>(n);
+        xt::xarray<float> y = A * xt::sin(2 * xt::numeric_constants<float>::PI * x * kf / nf);
+        auto res = xt::fft::fft(y) / static_cast<float>(n / 2);
         REQUIRE(A == doctest::Approx(std::abs(res(k))).epsilon(.0001));
     }
 
     TEST(xfft, ifft_power_2)
     {
         size_t k = 2;
         size_t n = 8;
-        size_t A = 10;
+        float A = 10;
         auto x = xt::linspace<float>(0, static_cast<float>(n - 1), n);
-        xt::xarray<float> y = A * xt::sin(2 * xt::numeric_constants<float>::PI * x * k / n);
-        auto res = xt::fft::ifft(y) / (n / 2);
+        const auto kf = static_cast<float>(k);
+        const auto nf = static_cast<float>(n);
+        xt::xarray<float> y = A * xt::sin(2 * xt::numeric_constants<float>::PI * x * kf / nf);
+        auto res = xt::fft::ifft(y) / static_cast<float>(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<float>(0, n - 1, n) * xt::ones<float>({dim, n});
-        xt::xarray<float> y = A * xt::sin(2 * xt::numeric_constants<float>::PI * x * k / n);
+        auto x = xt::linspace<float>(0, static_cast<float>(n - 1), n) * xt::ones<float>({dim, n});
+        const auto kf = static_cast<float>(k);
+        const auto nf = static_cast<float>(n);
+        xt::xarray<float> y = A * xt::sin(2 * xt::numeric_constants<float>::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<double>(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<float>(0, n - 1, n) * xt::ones<float>({dim, n});
-        xt::xarray<float> y = A * xt::sin(2 * xt::numeric_constants<float>::PI * x * k / n);
-        auto res = xt::fft::fft(y) / (n / 2.0);
+        auto x = xt::linspace<float>(0, static_cast<float>(n - 1), n) * xt::ones<float>({dim, n});
+        const auto kf = static_cast<float>(k);
+        const auto nf = static_cast<float>(n);
+        xt::xarray<float> y = A * xt::sin(2 * xt::numeric_constants<float>::PI * x * kf / nf);
+        auto res = xt::fft::fft(y) / (static_cast<double>(n) / 2);
         REQUIRE(A == doctest::Approx(std::abs(res(0, k))).epsilon(.0001));
         REQUIRE(A == doctest::Approx(std::abs(res(1, k))).epsilon(.0001));
     }