Repository: xtensor-stack/xtensor (language: python)

Issue:
filter missing broadcasting? view missing selector?
My eye crossed https://stackoverflow.com/q/74100940/2646505 . The goal there is to mimic
```python
import numpy as np

keep = np.array([True, False, True])
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
print(A[keep, :])
```

In xtensor this is not possible it seems.

## `xt::filter` does not broadcast:
```cpp
#include <xtensor/xarray.hpp>
#include <xtensor/xtensor.hpp>
#include <xtensor/xsort.hpp>
#include <xtensor/xview.hpp>
#include <xtensor/xindex_view.hpp>
#include <xtensor/xio.hpp>

int main()
{
    xt::xtensor<bool, 1> keep = {{true, false, true}};
    xt::xtensor<int, 2> A = {{1, 2, 3}, {4, 5, 6}, {4, 5, 6}};
    auto c = xt::filter(A, keep);
    std::cout << c << std::endl;
    return 0;
}
```

Returns only the first column!?
```
{1, 3}
```

## `xt::view` does not take logical

A syntax could be
```
int main()
{
    xt::xtensor<bool, 1> keep = {{true, false, true}};
    xt::xtensor<int, 2> A = {{1, 2, 3}, {4, 5, 6}, {4, 5, 6}};
    auto c = xt::view(A, xt::is_true(keep));
    std::cout << c << std::endl;
    return 0;
}
```
or have an overload on `keep`, but that might be a bit risky. What is the input is meant as logical, but given as `int`s?

## (bonus) View should take `xt::missing`

I realised that it would be nice to do
```python
import numpy as np

keep = np.array([True, False, True])
A = np.arange(3 * 4 * 5).reshape(3, 4, 5)
print(A[keep, ...])
```

as follows:
```cpp
int main()
{
    xt::xtensor<bool, 1> keep = {{true, false, true}};
    xt::xtensor<int, 2> A = {{1, 2, 3}, {4, 5, 6}, {4, 5, 6}};
    auto c = xt::view(A, xt::is_true(keep), xt::missing);
    std::cout << c << std::endl;
    return 0;
}
```



Discussion / hints:
Yes, `filter` doesn't take broadcasting. I 've talked about the implementation of `filter` in https://github.com/xtensor-stack/xtensor/issues/2584, and IMO, it's natural to have broadcasting in `filter`. I 'm curious about the intent of this implementation.
Sorry, I missed https://github.com/xtensor-stack/xtensor/issues/2584
Hi,

I think it would be nice to have this over any possible axis, as in NumPy
```cpp
auto c = xt::view(A, xt::drop(1), xt::keep(keep), xt::all());
```

Where `keep` is a dynamic array of bool the size of the relevant dimension, or a subset of indices to keep.
Would broadcasting be able to solve the general (boolean) case?

Fix patch:
diff --git a/docs/source/numpy.rst b/docs/source/numpy.rst
index ac2b0bebb..897d0383a 100644
--- a/docs/source/numpy.rst
+++ b/docs/source/numpy.rst
@@ -399,34 +399,40 @@ Comparisons
 Minimum, Maximum, Sorting
 -------------------------
 
-.. table::
+.. list-table::
    :widths: 50 50
+   :header-rows: 1
 
-   +-----------------------------------------------------+---------------------------------------------------------+
-   |            Python 3 - NumPy                         |                C++ 14 - xtensor                         |
-   +=====================================================+=========================================================+
-   | :any:`np.amin(a) <numpy.amin>`                      | :cpp:func:`xt::amin(a) <xt::amin>`                      |
-   +-----------------------------------------------------+---------------------------------------------------------+
-   | :any:`np.amax(a) <numpy.amax>`                      | :cpp:func:`xt::amax(a) <xt::amax>`                      |
-   +-----------------------------------------------------+---------------------------------------------------------+
-   | :any:`np.argmin(a) <numpy.argmin>`                  | :cpp:func:`xt::argmin(a) <xt::argmin>`                  |
-   +-----------------------------------------------------+---------------------------------------------------------+
-   | :any:`np.argmax(a, axis=1) <numpy.argmax>`          | :cpp:func:`xt::argmax(a, 1) <xt::argmax>`               |
-   +-----------------------------------------------------+---------------------------------------------------------+
-   | :any:`np.sort(a, axis=1) <numpy.sort>`              | :cpp:func:`xt::sort(a, 1) <xt::sort>`                   |
-   +-----------------------------------------------------+---------------------------------------------------------+
-   | :any:`np.argsort(a, axis=1) <numpy.argsort>`        | :cpp:func:`xt::argsort(a, 1) <xt::argsort>`             |
-   +-----------------------------------------------------+---------------------------------------------------------+
-   | :any:`np.unique(a) <numpy.unique>`                  | :cpp:func:`xt::unique(a) <xt::unique>`                  |
-   +-----------------------------------------------------+---------------------------------------------------------+
-   | :any:`np.setdiff1d(ar1, ar2) <numpy.setdiff1d>`     | :cpp:func:`xt::setdiff1d(ar1, ar2) <xt::setdiff1d>`     |
-   +-----------------------------------------------------+---------------------------------------------------------+
-   | :any:`np.partition(a, kth) <numpy.partition>`       | :cpp:func:`xt::partition(a, kth) <xt::partition>`       |
-   +-----------------------------------------------------+---------------------------------------------------------+
-   | :any:`np.argpartition(a, kth) <numpy.argpartition>` | :cpp:func:`xt::argpartition(a, kth) <xt::argpartition>` |
-   +-----------------------------------------------------+---------------------------------------------------------+
-   | :any:`np.median(a, axis) <numpy.median>`            | :cpp:func:`xt::median(a, axis) <xt::median>`            |
-   +-----------------------------------------------------+---------------------------------------------------------+
+   * - Python3 - NumPy
+     - C++14 - xtensor
+   * - :any:`np.amin(a) <numpy.amin>`
+     - :cpp:func:`xt::amin(a) <xt::amin>`
+   * - :any:`np.amax(a) <numpy.amax>`
+     - :cpp:func:`xt::amax(a) <xt::amax>`
+   * - :any:`np.argmin(a) <numpy.argmin>`
+     - :cpp:func:`xt::argmin(a) <xt::argmin>`
+   * - :any:`np.argmax(a, axis=1) <numpy.argmax>`
+     - :cpp:func:`xt::argmax(a, 1) <xt::argmax>`
+   * - :any:`np.sort(a, axis=1) <numpy.sort>`
+     - :cpp:func:`xt::sort(a, 1) <xt::sort>`
+   * - :any:`np.argsort(a, axis=1) <numpy.argsort>`
+     - :cpp:func:`xt::argsort(a, 1) <xt::argsort>`
+   * - :any:`np.unique(a) <numpy.unique>`
+     - :cpp:func:`xt::unique(a) <xt::unique>`
+   * - :any:`np.setdiff1d(ar1, ar2) <numpy.setdiff1d>`
+     - :cpp:func:`xt::setdiff1d(ar1, ar2) <xt::setdiff1d>`
+   * - :any:`np.partition(a, kth) <numpy.partition>`
+     - :cpp:func:`xt::partition(a, kth) <xt::partition>`
+   * - :any:`np.argpartition(a, kth) <numpy.argpartition>`
+     - :cpp:func:`xt::argpartition(a, kth) <xt::argpartition>`
+   * - :any:`np.quantile(a, [.1 .3], method="linear") <numpy.quantile>`
+     - :cpp:func:`xt::quantile(a, {.1, .3}, xt::quantile_method::linear) <xt::quantile>`
+   * - :any:`np.quantile(a, [.1, .3], axis=1 method="linear") <numpy.quantile>`
+     - :cpp:func:`xt::quantile(a, {.1, .3}, 1, xt::quantile_method::linear) <xt::quantile>`
+   * -
+     - :cpp:func:`xt::quantile(a, {.1, .3}, 1, 1.0, 1.0) <xt::quantile>`
+   * - :any:`np.median(a, axis=1) <numpy.median>`
+     - :cpp:func:`xt::median(a, 1) <xt::median>`
 
 Complex numbers
 ---------------
diff --git a/include/xtensor/xiterator.hpp b/include/xtensor/xiterator.hpp
index 55ee46ce9..afd77f73c 100644
--- a/include/xtensor/xiterator.hpp
+++ b/include/xtensor/xiterator.hpp
@@ -17,6 +17,7 @@
 #include <numeric>
 #include <vector>
 
+#include <xtl/xcompare.hpp>
 #include <xtl/xiterator_base.hpp>
 #include <xtl/xmeta_utils.hpp>
 #include <xtl/xsequence.hpp>
@@ -645,7 +646,7 @@ namespace xt
         {
             --i;
             size_type inc = (i == leading_i) ? n : 1;
-            if (index[i] + inc < shape[i])
+            if (xtl::cmp_less(index[i] + inc, shape[i]))
             {
                 index[i] += inc;
                 stepper.step(i, inc);
@@ -723,7 +724,7 @@ namespace xt
         {
             --i;
             size_type inc = (i == leading_i) ? n : 1;
-            if (index[i] >= inc)
+            if (xtl::cmp_greater_equal(index[i], inc))
             {
                 index[i] -= inc;
                 stepper.step_back(i, inc);
diff --git a/include/xtensor/xsort.hpp b/include/xtensor/xsort.hpp
index 389a3c1fd..a8253dd55 100644
--- a/include/xtensor/xsort.hpp
+++ b/include/xtensor/xsort.hpp
@@ -11,14 +11,22 @@
 #define XTENSOR_SORT_HPP
 
 #include <algorithm>
+#include <cmath>
 #include <utility>
 
+#include <xtl/xcompare.hpp>
+
+#include "xadapt.hpp"
 #include "xarray.hpp"
 #include "xeval.hpp"
+#include "xindex_view.hpp"
 #include "xmanipulation.hpp"
+#include "xmath.hpp"
 #include "xslice.hpp"  // for xnone
 #include "xtensor.hpp"
 #include "xtensor_config.hpp"
+#include "xtensor_forward.hpp"
+#include "xview.hpp"
 
 namespace xt
 {
@@ -772,6 +780,328 @@ namespace xt
         return argpartition(e, std::array<std::size_t, 1>({kth}), axis);
     }
 
+    /******************
+     *  xt::quantile  *
+     ******************/
+
+    namespace detail
+    {
+        template <class S, class I, class K, class O>
+        inline void select_indices_impl(
+            const S& shape,
+            const I& indices,
+            std::size_t axis,
+            std::size_t current_dim,
+            const K& current_index,
+            O& out
+        )
+        {
+            using id_t = typename K::value_type;
+            if ((current_dim < shape.size() - 1) && (current_dim == axis))
+            {
+                for (auto i : indices)
+                {
+                    auto idx = current_index;
+                    idx[current_dim] = i;
+                    select_indices_impl(shape, indices, axis, current_dim + 1, idx, out);
+                }
+            }
+            else if ((current_dim < shape.size() - 1) && (current_dim != axis))
+            {
+                for (id_t i = 0; xtl::cmp_less(i, shape[current_dim]); ++i)
+                {
+                    auto idx = current_index;
+                    idx[current_dim] = i;
+                    select_indices_impl(shape, indices, axis, current_dim + 1, idx, out);
+                }
+            }
+            else if ((current_dim == shape.size() - 1) && (current_dim == axis))
+            {
+                for (auto i : indices)
+                {
+                    auto idx = current_index;
+                    idx[current_dim] = i;
+                    out.push_back(std::move(idx));
+                }
+            }
+            else if ((current_dim == shape.size() - 1) && (current_dim != axis))
+            {
+                for (id_t i = 0; xtl::cmp_less(i, shape[current_dim]); ++i)
+                {
+                    auto idx = current_index;
+                    idx[current_dim] = i;
+                    out.push_back(std::move(idx));
+                }
+            }
+        }
+
+        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>;
+            auto out = std::vector<index_type>();
+            select_indices_impl(shape, indices, axis, 0, xtl::make_sequence<index_type>(shape.size()), out);
+            return out;
+        }
+
+        // TODO remove when fancy index views are implemented
+        // Poor man's indexing along a single axis as in NumPy a[:, [1, 3, 4]]
+        template <class E, class I>
+        inline auto fancy_indexing(E&& e, const I& indices, std::ptrdiff_t axis)
+        {
+            std::size_t const ax = normalize_axis(e.dimension(), axis);
+            using shape_t = get_strides_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(
+                index_view(std::forward<E>(e), select_indices(e.shape(), indices, ax)),
+                std::move(shape)
+            );
+        }
+
+        template <class T, class I, class P>
+        inline auto quantile_kth_gamma(std::size_t n, const P& probas, T alpha, T beta)
+        {
+            const auto m = alpha + probas * (T(1) - alpha - beta);
+            // Evaluting since reused a lot
+            const auto p_n_m = eval(probas * static_cast<T>(n) + m - 1);
+            // Previous (virtual) index, may be out of bounds
+            const auto j = floor(p_n_m);
+            const auto j_jp1 = concatenate(xtuple(j, j + 1));
+            // Both interpolation indices, k and k+1
+            const auto k_kp1 = xt::cast<std::size_t>(clip(j_jp1, T(0), T(n - 1)));
+            // Both interpolation coefficients, 1-gamma and gamma
+            const auto omg_g = concatenate(xtuple(T(1) - (p_n_m - j), p_n_m - j));
+            return std::make_pair(eval(k_kp1), eval(omg_g));
+        }
+
+        // TODO should implement unsqueeze rather
+        template <class S>
+        inline auto unsqueeze_shape(const S& shape, std::size_t axis)
+        {
+            XTENSOR_ASSERT(axis <= shape.size());
+            auto new_shape = xtl::forward_sequence<xt::svector<std::size_t>, decltype(shape)>(shape);
+            new_shape.insert(new_shape.begin() + axis, 1);
+            return new_shape;
+        }
+    }
+
+    /**
+     * Compute quantiles over the given axis.
+     *
+     * In a sorted array represneting a distribution of numbers, the quantile of a probability ``p``
+     * is the the cut value ``q`` such that a fraction ``p`` of the distribution is lesser or equal
+     * to ``q``.
+     * When the cutpoint falls between two elemnts of the sample distribution, a interpolation is
+     * computed using the @p alpha and @p beta coefficients, as descripted in
+     * (Hyndman and Fan, 1996).
+     *
+     * The algorithm partially sorts entries in a copy along the @p axis axis.
+     *
+     * @ingroup xt_xsort
+     * @param e Expression containing the distribution over which the quantiles are computed.
+     * @param probas An list of probability associated with each desired quantiles.
+     *        All elements must be in the range ``[0, 1]``.
+     * @param axis The dimension in which to compute the quantiles, *i.e* the axis representing the
+     *        distribution.
+     * @param alpha Interpolation parameter. Must be in the range ``[0, 1]]``.
+     * @param beta Interpolation parameter. Must be in the range ``[0, 1]]``.
+     * @tparam T The type in which the quatile are computed.
+     * @return An expression with as many dimensions as the input @p e.
+     *         The first axis correspond to the quantiles.
+     *         The other axes are the axes that remain after the reduction of @p e.
+     * @see (Hyndman and Fan, 1996) R. J. Hyndman and Y. Fan,
+     *      "Sample quantiles in statistical packages", The American Statistician,
+     *      50(4), pp. 361-365, 1996
+     * @see https://en.wikipedia.org/wiki/Quantile
+     */
+    template <class T = double, class E, class P>
+    inline auto quantile(E&& e, const P& probas, std::ptrdiff_t axis, T alpha, T beta)
+    {
+        XTENSOR_ASSERT(all(0. <= probas));
+        XTENSOR_ASSERT(all(probas <= 1.));
+        XTENSOR_ASSERT(0. <= alpha);
+        XTENSOR_ASSERT(alpha <= 1.);
+        XTENSOR_ASSERT(0. <= beta);
+        XTENSOR_ASSERT(beta <= 1.);
+
+        using tmp_shape_t = get_strides_t<typename std::decay_t<E>::shape_type>;
+        using id_t = typename tmp_shape_t::value_type;
+
+        std::size_t const ax = normalize_axis(e.dimension(), axis);
+        std::size_t const 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);
+
+        // Reshape interpolation coefficients
+        auto gm1_g_shape = xtl::make_sequence<tmp_shape_t>(e.dimension(), 1);
+        gm1_g_shape[ax] = kth_gamma.second.size();
+        auto gm1_g_reshaped = reshape_view(std::move(kth_gamma.second), std::move(gm1_g_shape));
+
+        // Compute interpolation
+        // TODO(C++20) use (and create) xt::lerp in C++
+        auto e_kth_g = std::move(e_kth) * std::move(gm1_g_reshaped);
+        // Reshape pairwise interpolate for suming along new axis
+        auto e_kth_g_shape = detail::unsqueeze_shape(e_kth_g.shape(), ax);
+        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);
+    }
+
+    // Static proba array overload
+    template <class T = double, class E, std::size_t N>
+    inline auto quantile(E&& e, const T (&probas)[N], std::ptrdiff_t axis, T alpha, T beta)
+    {
+        return quantile(std::forward<E>(e), adapt(probas, {N}), axis, alpha, beta);
+    }
+
+    /**
+     * Compute quantiles of the whole expression.
+     *
+     * The quantiles are computed over the whole expression, as if flatten in a one-dimensional
+     * expression.
+     *
+     * @ingroup xt_xsort
+     * @see xt::quantile(E&& e, P const& probas, std::ptrdiff_t axis, T alpha, T beta)
+     */
+    template <class T = double, class E, class P>
+    inline auto quantile(E&& e, const P& probas, T alpha, T beta)
+    {
+        return quantile(xt::ravel(std::forward<E>(e)), probas, 0, alpha, beta);
+    }
+
+    // Static proba array overload
+    template <class T = double, class E, std::size_t N>
+    inline auto quantile(E&& e, const T (&probas)[N], T alpha, T beta)
+    {
+        return quantile(std::forward<E>(e), adapt(probas, {N}), alpha, beta);
+    }
+
+    /**
+     * Quantile interpolation method.
+     *
+     * Predefined methods for interpolating quantiles, as defined in (Hyndman and Fan, 1996).
+     *
+     * @ingroup xt_xsort
+     * @see (Hyndman and Fan, 1996) R. J. Hyndman and Y. Fan,
+     *      "Sample quantiles in statistical packages", The American Statistician,
+     *      50(4), pp. 361-365, 1996
+     * @see xt::quantile(E&& e, P const& probas, std::ptrdiff_t axis, xt::quantile_method method)
+     */
+    enum class quantile_method
+    {
+        /** Method 4 of (Hyndman and Fan, 1996) with ``alpha=0`` and ``beta=1``. */
+        interpolated_inverted_cdf = 4,
+        /** Method 5 of (Hyndman and Fan, 1996) with ``alpha=1/2`` and ``beta=1/2``. */
+        hazen,
+        /** Method 6 of (Hyndman and Fan, 1996) with ``alpha=0`` and ``beta=0``. */
+        weibull,
+        /** Method 7 of (Hyndman and Fan, 1996) with ``alpha=1`` and ``beta=1``. */
+        linear,
+        /** Method 8 of (Hyndman and Fan, 1996) with ``alpha=1/3`` and ``beta=1/3``. */
+        median_unbiased,
+        /** Method 9 of (Hyndman and Fan, 1996) with ``alpha=3/8`` and ``beta=3/8``. */
+        normal_unbiased,
+    };
+
+    /**
+     * Compute quantiles over the given axis.
+     *
+     * The function takes the name of a predefined method to compute to interpolate between values.
+     *
+     * @ingroup xt_xsort
+     * @see xt::quantile_method
+     * @see xt::quantile(E&& e, P const& probas, std::ptrdiff_t axis, T alpha, T beta)
+     */
+    template <class T = double, class E, class P>
+    inline auto
+    quantile(E&& e, const P& probas, std::ptrdiff_t axis, quantile_method method = quantile_method::linear)
+    {
+        T alpha = 0.;
+        T beta = 0.;
+        switch (method)
+        {
+            case (quantile_method::interpolated_inverted_cdf):
+            {
+                alpha = 0.;
+                beta = 1.;
+                break;
+            }
+            case (quantile_method::hazen):
+            {
+                alpha = 0.5;
+                beta = 0.5;
+                break;
+            }
+            case (quantile_method::weibull):
+            {
+                alpha = 0.;
+                beta = 0.;
+                break;
+            }
+            case (quantile_method::linear):
+            {
+                alpha = 1.;
+                beta = 1.;
+                break;
+            }
+            case (quantile_method::median_unbiased):
+            {
+                alpha = 1. / 3.;
+                beta = 1. / 3.;
+                break;
+            }
+            case (quantile_method::normal_unbiased):
+            {
+                alpha = 3. / 8.;
+                beta = 3. / 8.;
+                break;
+            }
+        }
+        return quantile(std::forward<E>(e), probas, axis, alpha, beta);
+    }
+
+    // Static proba array overload
+    template <class T = double, class E, std::size_t N>
+    inline auto
+    quantile(E&& e, const T (&probas)[N], std::ptrdiff_t axis, quantile_method method = quantile_method::linear)
+    {
+        return quantile(std::forward<E>(e), adapt(probas, {N}), axis, method);
+    }
+
+    /**
+     * Compute quantiles of the whole expression.
+     *
+     * The quantiles are computed over the whole expression, as if flatten in a one-dimensional
+     * expression.
+     * The function takes the name of a predefined method to compute to interpolate between values.
+     *
+     * @ingroup xt_xsort
+     * @see xt::quantile_method
+     * @see xt::quantile(E&& e, P const& probas, std::ptrdiff_t axis, xt::quantile_method method)
+     */
+    template <class T = double, class E, class P>
+    inline auto quantile(E&& e, const P& probas, quantile_method method = quantile_method::linear)
+    {
+        return quantile(xt::ravel(std::forward<E>(e)), probas, 0, method);
+    }
+
+    // Static proba array overload
+    template <class T = double, class E, std::size_t N>
+    inline auto quantile(E&& e, const T (&probas)[N], quantile_method method = quantile_method::linear)
+    {
+        return quantile(std::forward<E>(e), adapt(probas, {N}), method);
+    }
+
+    /****************
+     *  xt::median  *
+     ****************/
+
     template <class E>
     inline typename std::decay_t<E>::value_type median(E&& e)
     {


Test patch:
diff --git a/test/files/cppy_source/test_extended_xsort.cppy b/test/files/cppy_source/test_extended_xsort.cppy
index 3d2581f96..ab89affe9 100644
--- a/test/files/cppy_source/test_extended_xsort.cppy
+++ b/test/files/cppy_source/test_extended_xsort.cppy
@@ -11,6 +11,7 @@
 #include <xtensor/xio.hpp>
 #include <xtensor/xview.hpp>
 #include <xtensor/xsort.hpp>
+#include <xtensor/xmath.hpp>
 
 #include "test_common_macros.hpp"
 
@@ -128,6 +129,65 @@ namespace xt
         EXPECT_EQ(py_a[part_a1(12)], py_a[static_cast<std::size_t>(py_a1(12))]);
     }
 
+    /*py
+    a = np.random.rand(4, 5, 6)
+    */
+    TEST(xtest_extended, quantile_axis)
+    {
+        // py_a
+
+        // py_q0 = np.quantile(a, [0., .3, .1, 1.], axis=0)
+        // py_q1 = np.quantile(a, [0., .3, .1, 1.], axis=1)
+        // py_q2 = np.quantile(a, [0., .3, .1, 1.], axis=2)
+
+        EXPECT_TRUE(xt::allclose(py_q0, xt::quantile(py_a, {0., .3, .1, 1.}, 0)));
+        EXPECT_TRUE(xt::allclose(py_q1, xt::quantile(py_a, {0., .3, .1, 1.}, 1)));
+        EXPECT_TRUE(xt::allclose(py_q2, xt::quantile(py_a, {0., .3, .1, 1.}, 2)));
+    }
+
+    /*py
+    a = np.random.rand(4, 5, 6)
+    */
+    TEST(xtest_extended, quantile_methods)
+    {
+        // py_a
+
+        // py_q4 = np.quantile(a, [0., .3, .1, 1.], method="interpolated_inverted_cdf")
+        // py_q5 = np.quantile(a, [0., .3, .1, 1.], method="hazen")
+        // py_q6 = np.quantile(a, [0., .3, .1, 1.], method="weibull")
+        // py_q7 = np.quantile(a, [0., .3, .1, 1.], method="linear")
+        // py_q8 = np.quantile(a, [0., .3, .1, 1.], method="median_unbiased")
+        // py_q9 = np.quantile(a, [0., .3, .1, 1.], method="normal_unbiased")
+
+        EXPECT_TRUE(
+            xt::allclose(
+                py_q4,
+                xt::quantile(py_a, {0., .3, .1, 1.}, quantile_method::interpolated_inverted_cdf)
+            )
+        );
+        EXPECT_TRUE(
+            xt::allclose(py_q5, xt::quantile(py_a, {0., .3, .1, 1.}, quantile_method::hazen))
+        );
+        EXPECT_TRUE(
+            xt::allclose(py_q6, xt::quantile(py_a, {0., .3, .1, 1.}, quantile_method::weibull))
+        );
+        EXPECT_TRUE(
+            xt::allclose(py_q7, xt::quantile(py_a, {0., .3, .1, 1.}, quantile_method::linear))
+        );
+        EXPECT_TRUE(
+            xt::allclose(
+                py_q8,
+                xt::quantile(py_a, {0., .3, .1, 1.}, quantile_method::median_unbiased)
+            )
+        );
+        EXPECT_TRUE(
+            xt::allclose(
+                py_q9,
+                xt::quantile(py_a, {0., .3, .1, 1.}, quantile_method::normal_unbiased)
+            )
+        );
+    }
+
     /*py
     a = np.random.rand(5, 5, 5)
     */
diff --git a/test/test_extended_xsort.cpp b/test/test_extended_xsort.cpp
index fd75604a4..b17fb8837 100644
--- a/test/test_extended_xsort.cpp
+++ b/test/test_extended_xsort.cpp
@@ -12,6 +12,7 @@
 
 #include <xtensor/xarray.hpp>
 #include <xtensor/xio.hpp>
+#include <xtensor/xmath.hpp>
 #include <xtensor/xsort.hpp>
 #include <xtensor/xview.hpp>
 
@@ -107,7 +108,7 @@ namespace xt
 
         // py_median = np.median(a)
         double py_median = 300.0;
-        EXPECT_EQ(static_cast<std::decay<decltype(py_a)>::type::value_type>(py_median), xt::median(py_a));
+        EXPECT_EQ(static_cast<decltype(py_a)::value_type>(py_median), xt::median(py_a));
     }
 
     /*py
@@ -151,6 +152,544 @@ namespace xt
         EXPECT_EQ(py_a[part_a1(12)], py_a[static_cast<std::size_t>(py_a1(12))]);
     }
 
+    /*py
+    a = np.random.rand(4, 5, 6)
+    */
+    TEST(xtest_extended, quantile_axis)
+    {
+        // py_a
+        xarray<double> py_a = {
+            {{0.0650515929852795,
+              0.9488855372533332,
+              0.9656320330745594,
+              0.8083973481164611,
+              0.3046137691733707,
+              0.0976721140063839},
+             {0.6842330265121569,
+              0.4401524937396013,
+              0.1220382348447788,
+              0.4951769101112702,
+              0.0343885211152184,
+              0.9093204020787821},
+             {0.2587799816000169,
+              0.662522284353982,
+              0.311711076089411,
+              0.5200680211778108,
+              0.5467102793432796,
+              0.184854455525527},
+             {0.9695846277645586,
+              0.7751328233611146,
+              0.9394989415641891,
+              0.8948273504276488,
+              0.5978999788110851,
+              0.9218742350231168},
+             {0.0884925020519195,
+              0.1959828624191452,
+              0.0452272889105381,
+              0.3253303307632643,
+              0.388677289689482,
+              0.2713490317738959}},
+
+            {{0.8287375091519293,
+              0.3567533266935893,
+              0.2809345096873808,
+              0.5426960831582485,
+              0.1409242249747626,
+              0.8021969807540397},
+             {0.0745506436797708,
+              0.9868869366005173,
+              0.7722447692966574,
+              0.1987156815341724,
+              0.0055221171236024,
+              0.8154614284548342},
+             {0.7068573438476171,
+              0.7290071680409873,
+              0.7712703466859457,
+              0.0740446517340904,
+              0.3584657285442726,
+              0.1158690595251297},
+             {0.8631034258755935,
+              0.6232981268275579,
+              0.3308980248526492,
+              0.0635583502860236,
+              0.3109823217156622,
+              0.325183322026747},
+             {0.7296061783380641,
+              0.6375574713552131,
+              0.8872127425763265,
+              0.4722149251619493,
+              0.1195942459383017,
+              0.713244787222995}},
+
+            {{0.7607850486168974,
+              0.5612771975694962,
+              0.770967179954561,
+              0.4937955963643907,
+              0.5227328293819941,
+              0.4275410183585496},
+             {0.0254191267440952,
+              0.1078914269933045,
+              0.0314291856867343,
+              0.6364104112637804,
+              0.3143559810763267,
+              0.5085706911647028},
+             {0.907566473926093,
+              0.2492922291488749,
+              0.4103829230356297,
+              0.7555511385430487,
+              0.2287981654916225,
+              0.076979909828793},
+             {0.289751452913768,
+              0.1612212872540044,
+              0.9296976523425731,
+              0.808120379564417,
+              0.6334037565104235,
+              0.8714605901877177},
+             {0.8036720768991145,
+              0.1865700588860358,
+              0.8925589984899778,
+              0.5393422419156507,
+              0.8074401551640625,
+              0.8960912999234932}},
+
+            {{0.3180034749718639,
+              0.1100519245276768,
+              0.2279351625419417,
+              0.4271077886262563,
+              0.8180147659224931,
+              0.8607305832563434},
+             {0.0069521305311907,
+              0.5107473025775657,
+              0.417411003148779,
+              0.2221078104707302,
+              0.1198653673336828,
+              0.337615171403628},
+             {0.9429097039125192,
+              0.3232029320207552,
+              0.5187906217433661,
+              0.7030189588951778,
+              0.363629602379294,
+              0.9717820827209607},
+             {0.9624472949421112,
+              0.2517822958253642,
+              0.4972485058923855,
+              0.3008783098167697,
+              0.2848404943774676,
+              0.0368869473545328},
+             {0.6095643339798968,
+              0.5026790232288615,
+              0.0514787512499894,
+              0.2786464642366114,
+              0.9082658859666537,
+              0.2395618906669724}}};
+
+        // py_q0 = np.quantile(a, [0., .3, .1, 1.], axis=0)
+        xarray<double> py_q0 = {
+            {{0.0650515929852795,
+              0.1100519245276768,
+              0.2279351625419417,
+              0.4271077886262563,
+              0.1409242249747626,
+              0.0976721140063839},
+             {0.0069521305311907,
+              0.1078914269933045,
+              0.0314291856867343,
+              0.1987156815341724,
+              0.0055221171236024,
+              0.337615171403628},
+             {0.2587799816000169,
+              0.2492922291488749,
+              0.311711076089411,
+              0.0740446517340904,
+              0.2287981654916225,
+              0.076979909828793},
+             {0.289751452913768,
+              0.1612212872540044,
+              0.3308980248526492,
+              0.0635583502860236,
+              0.2848404943774676,
+              0.0368869473545328},
+             {0.0884925020519195,
+              0.1865700588860358,
+              0.0452272889105381,
+              0.2786464642366114,
+              0.1195942459383017,
+              0.2395618906669724}},
+
+            {{0.2927082867732054,
+              0.332083186476998,
+              0.2756345749728368,
+              0.4871268155905773,
+              0.2882448147535099,
+              0.394554127923333},
+             {0.0235724271228047,
+              0.4069263870649716,
+              0.1129773299289744,
+              0.2197685975770745,
+              0.0315018807160568,
+              0.4914751391885953},
+             {0.662049607622857,
+              0.3158118617335672,
+              0.4005157383410078,
+              0.4754656842334387,
+              0.3454989722390076,
+              0.111980144555496},
+             {0.8057682285794109,
+              0.2427261949682282,
+              0.4806134577884118,
+              0.277146313863695,
+              0.3083681389818427,
+              0.2963536845595256},
+             {0.5574571507870991,
+              0.1950415820658343,
+              0.0508536050160442,
+              0.3206619441105991,
+              0.361768985314364,
+              0.2681703176632035}},
+
+            {{0.1409371575812548,
+              0.1840623451774505,
+              0.2438349666855734,
+              0.4471141309476966,
+              0.1900310882343451,
+              0.1966327853120336},
+             {0.012492229395062,
+              0.2075697470171935,
+              0.0586119004341476,
+              0.2057333202151398,
+              0.0141820383210872,
+              0.3889018273319504},
+             {0.393203190274297,
+              0.271465440010439,
+              0.3413126301732766,
+              0.2078516625672065,
+              0.2676984344074175,
+              0.088646654737694},
+             {0.4617570448023157,
+              0.1883895898254123,
+              0.3808031691645701,
+              0.1347543381452475,
+              0.292683042578926,
+              0.1233758597561971},
+             {0.2448140516303127,
+              0.1893938999459686,
+              0.0471027276123735,
+              0.2926516241946073,
+              0.2003191590636558,
+              0.2490980329990495}},
+
+            {{0.8287375091519293,
+              0.9488855372533332,
+              0.9656320330745594,
+              0.8083973481164611,
+              0.8180147659224931,
+              0.8607305832563434},
+             {0.6842330265121569,
+              0.9868869366005173,
+              0.7722447692966574,
+              0.6364104112637804,
+              0.3143559810763267,
+              0.9093204020787821},
+             {0.9429097039125192,
+              0.7290071680409873,
+              0.7712703466859457,
+              0.7555511385430487,
+              0.5467102793432796,
+              0.9717820827209607},
+             {0.9695846277645586,
+              0.7751328233611146,
+              0.9394989415641891,
+              0.8948273504276488,
+              0.6334037565104235,
+              0.9218742350231168},
+             {0.8036720768991145,
+              0.6375574713552131,
+              0.8925589984899778,
+              0.5393422419156507,
+              0.9082658859666537,
+              0.8960912999234932}}};
+        // py_q1 = np.quantile(a, [0., .3, .1, 1.], axis=1)
+        xarray<double> py_q1 = {
+            {{0.0650515929852795,
+              0.1959828624191452,
+              0.0452272889105381,
+              0.3253303307632643,
+              0.0343885211152184,
+              0.0976721140063839},
+             {0.0745506436797708,
+              0.3567533266935893,
+              0.2809345096873808,
+              0.0635583502860236,
+              0.0055221171236024,
+              0.1158690595251297},
+             {0.0254191267440952,
+              0.1078914269933045,
+              0.0314291856867343,
+              0.4937955963643907,
+              0.2287981654916225,
+              0.076979909828793},
+             {0.0069521305311907,
+              0.1100519245276768,
+              0.0514787512499894,
+              0.2221078104707302,
+              0.1198653673336828,
+              0.0368869473545328}},
+
+            {{0.122549997961539,
+              0.4846264518624774,
+              0.1599728030937052,
+              0.5001551323245783,
+              0.321426473276593,
+              0.2021533707752008},
+             {0.7114071107457065,
+              0.6261499957330889,
+              0.4189724892193085,
+              0.0989788576941068,
+              0.1238602417455939,
+              0.4027956150659966},
+             {0.3839581720543939,
+              0.1662910415804107,
+              0.4824997744194159,
+              0.5587558757852766,
+              0.3560313507374601,
+              0.4437469529197803},
+             {0.3763156467734705,
+              0.2660664230644424,
+              0.2658303306633091,
+              0.2830928333526431,
+              0.3005983159778329,
+              0.2591725468143035}},
+
+            {{0.0744279566119355,
+              0.2936507149473276,
+              0.0759516672842344,
+              0.3932689625024667,
+              0.1424786203384793,
+              0.1325450506140411},
+             {0.3274733237469094,
+              0.4633712467471767,
+              0.3009199157534881,
+              0.0677528708652503,
+              0.0511509686494821,
+              0.1995947645257767},
+             {0.1311520572119643,
+              0.1292233710975844,
+              0.1830106806262924,
+              0.5120142545848947,
+              0.2630212917255041,
+              0.2172043532406956},
+             {0.13137266830746,
+              0.1667440730467517,
+              0.1220613157667703,
+              0.2447232719770827,
+              0.1858554181511967,
+              0.1179569246795086}},
+
+            {{0.9695846277645586,
+              0.9488855372533332,
+              0.9656320330745594,
+              0.8948273504276488,
+              0.5978999788110851,
+              0.9218742350231168},
+             {0.8631034258755935,
+              0.9868869366005173,
+              0.8872127425763265,
+              0.5426960831582485,
+              0.3584657285442726,
+              0.8154614284548342},
+             {0.907566473926093,
+              0.5612771975694962,
+              0.9296976523425731,
+              0.808120379564417,
+              0.8074401551640625,
+              0.8960912999234932},
+             {0.9624472949421112,
+              0.5107473025775657,
+              0.5187906217433661,
+              0.7030189588951778,
+              0.9082658859666537,
+              0.9717820827209607}}};
+        // py_q2 = np.quantile(a, [0., .3, .1, 1.], axis=2)
+        xarray<double> py_q2 = {
+            {{0.0650515929852795, 0.0343885211152184, 0.184854455525527, 0.5978999788110851, 0.0452272889105381},
+             {0.1409242249747626, 0.0055221171236024, 0.0740446517340904, 0.0635583502860236, 0.1195942459383017},
+             {0.4275410183585496, 0.0254191267440952, 0.076979909828793, 0.1612212872540044, 0.1865700588860358},
+             {0.1100519245276768, 0.0069521305311907, 0.3232029320207552, 0.0368869473545328, 0.0514787512499894}},
+
+            {{0.2011429415898773, 0.2810953642921901, 0.2852455288447139, 0.8349800868943817, 0.1422376822355323},
+             {0.318843918190485, 0.1366331626069716, 0.2371673940347012, 0.3180828218712046, 0.5548861982585812},
+             {0.5082642128731925, 0.0696603063400194, 0.2390451973202487, 0.4615776047120957, 0.6715071594073826},
+             {0.2729693187569028, 0.1709865889022065, 0.44121011206133, 0.2683113951014159, 0.2591041774517919}},
+
+            {{0.0813618534958317, 0.0782133779799986, 0.221817218562772, 0.6865164010860998, 0.0668598954812288},
+             {0.2109293673310717, 0.0400363804016866, 0.09495685562961, 0.1872703360008429, 0.2959045855501255},
+             {0.4606683073614702, 0.0284241562154147, 0.1528890376602077, 0.2254863700838862, 0.3629561504008433},
+             {0.1689935435348092, 0.0634087489324368, 0.3434162672000246, 0.1443346215899485, 0.1455203209584809}},
+
+            {{0.9656320330745594, 0.9093204020787821, 0.662522284353982, 0.9695846277645586, 0.388677289689482},
+             {0.8287375091519293, 0.9868869366005173, 0.7712703466859457, 0.8631034258755935, 0.8872127425763265},
+             {0.770967179954561, 0.6364104112637804, 0.907566473926093, 0.9296976523425731, 0.8960912999234932},
+             {0.8607305832563434, 0.5107473025775657, 0.9717820827209607, 0.9624472949421112, 0.9082658859666537}}};
+
+        EXPECT_TRUE(xt::allclose(py_q0, xt::quantile(py_a, {0., .3, .1, 1.}, 0)));
+        EXPECT_TRUE(xt::allclose(py_q1, xt::quantile(py_a, {0., .3, .1, 1.}, 1)));
+        EXPECT_TRUE(xt::allclose(py_q2, xt::quantile(py_a, {0., .3, .1, 1.}, 2)));
+    }
+
+    /*py
+    a = np.random.rand(4, 5, 6)
+    */
+    TEST(xtest_extended, quantile_methods)
+    {
+        // py_a
+        xarray<double> py_a = {
+            {{0.1448948720912231,
+              0.489452760277563,
+              0.9856504541106007,
+              0.2420552715115004,
+              0.6721355474058786,
+              0.7616196153287176},
+             {0.2376375439923997,
+              0.7282163486118596,
+              0.3677831327192532,
+              0.6323058305935795,
+              0.6335297107608947,
+              0.5357746840747585},
+             {0.0902897700544083,
+              0.835302495589238,
+              0.3207800649717358,
+              0.1865185103998542,
+              0.0407751415547639,
+              0.5908929431882418},
+             {0.6775643618422824,
+              0.0165878289278562,
+              0.512093058299281,
+              0.226495775197938,
+              0.6451727904094499,
+              0.1743664290049914},
+             {0.690937738102466,
+              0.3867353463005374,
+              0.9367299887367345,
+              0.1375209441459933,
+              0.3410663510502585,
+              0.1134735212405891}},
+
+            {{0.9246936182785628,
+              0.877339353380981,
+              0.2579416277151556,
+              0.659984046034179,
+              0.8172222002012158,
+              0.5552008115994623},
+             {0.5296505783560065,
+              0.2418522909004517,
+              0.0931027678058992,
+              0.8972157579533268,
+              0.9004180571633305,
+              0.6331014572732679},
+             {0.3390297910487007,
+              0.3492095746126609,
+              0.7259556788702394,
+              0.8971102599525771,
+              0.8870864242651173,
+              0.7798755458576239},
+             {0.6420316461542878,
+              0.0841399649950488,
+              0.1616287140946138,
+              0.8985541885270792,
+              0.6064290596595899,
+              0.0091970516166296},
+             {0.1014715428660321,
+              0.6635017691080558,
+              0.0050615838462187,
+              0.1608080514174987,
+              0.5487337893665861,
+              0.6918951976926933}},
+
+            {{0.6519612595026005,
+              0.2242693094605598,
+              0.7121792213475359,
+              0.2372490874968001,
+              0.3253996981592677,
+              0.7464914051180241},
+             {0.6496328990472147,
+              0.8492234104941779,
+              0.6576128923003434,
+              0.5683086033354716,
+              0.0936747678280925,
+              0.3677158030594335},
+             {0.2652023676817254,
+              0.2439896433790836,
+              0.9730105547524456,
+              0.3930977246667604,
+              0.8920465551771133,
+              0.6311386259972629},
+             {0.7948113035416484,
+              0.5026370931051921,
+              0.5769038846263591,
+              0.4925176938188639,
+              0.1952429877980445,
+              0.7224521152615053},
+             {0.2807723624408558,
+              0.0243159664314538,
+              0.6454722959071678,
+              0.1771106794070489,
+              0.9404585843529143,
+              0.9539285770025874}},
+
+            {{0.9148643902204485,
+              0.3701587002554444,
+              0.0154566165288674,
+              0.9283185625877254,
+              0.4281841483173143,
+              0.9666548190436696},
+             {0.9636199770892528,
+              0.8530094554673601,
+              0.2944488920695857,
+              0.3850977286019253,
+              0.8511366715168569,
+              0.3169220051562777},
+             {0.1694927466860925,
+              0.5568012624583502,
+              0.936154774160781,
+              0.696029796674973,
+              0.570061170089365,
+              0.0971764937707685},
+             {0.6150072266991697,
+              0.9900538501042633,
+              0.140084015236524,
+              0.5183296523637367,
+              0.8773730719279554,
+              0.7407686177542044},
+             {0.697015740995268,
+              0.7024840839871093,
+              0.3594911512197552,
+              0.2935918442644934,
+              0.8093611554785136,
+              0.8101133946791808}}};
+
+        // py_q4 = np.quantile(a, [0., .3, .1, 1.], method="interpolated_inverted_cdf")
+        xarray<double> py_q4 = {0.0050615838462187, 0.3169220051562777, 0.1014715428660321, 0.9900538501042633};
+        // py_q5 = np.quantile(a, [0., .3, .1, 1.], method="hazen")
+        xarray<double> py_q5 = {0.0050615838462187, 0.3188510350640067, 0.1074725320533106, 0.9900538501042633};
+        // py_q6 = np.quantile(a, [0., .3, .1, 1.], method="weibull")
+        xarray<double> py_q6 = {0.0050615838462187, 0.3180794231009151, 0.1026717407034878, 0.9900538501042633};
+        // py_q7 = np.quantile(a, [0., .3, .1, 1.], method="linear")
+        xarray<double> py_q7 = {0.0050615838462187, 0.3196226470270984, 0.1122733234031334, 0.9900538501042633};
+        // py_q8 = np.quantile(a, [0., .3, .1, 1.], method="median_unbiased")
+        xarray<double> py_q8 = {0.0050615838462187, 0.3185938310763095, 0.1058722682700363, 0.9900538501042633};
+        // py_q9 = np.quantile(a, [0., .3, .1, 1.], method="normal_unbiased")
+        xarray<double> py_q9 = {0.0050615838462187, 0.3186581320732338, 0.1062723342158549, 0.9900538501042633};
+
+        EXPECT_TRUE(
+            xt::allclose(py_q4, xt::quantile(py_a, {0., .3, .1, 1.}, quantile_method::interpolated_inverted_cdf))
+        );
+        EXPECT_TRUE(xt::allclose(py_q5, xt::quantile(py_a, {0., .3, .1, 1.}, quantile_method::hazen)));
+        EXPECT_TRUE(xt::allclose(py_q6, xt::quantile(py_a, {0., .3, .1, 1.}, quantile_method::weibull)));
+        EXPECT_TRUE(xt::allclose(py_q7, xt::quantile(py_a, {0., .3, .1, 1.}, quantile_method::linear)));
+        EXPECT_TRUE(xt::allclose(py_q8, xt::quantile(py_a, {0., .3, .1, 1.}, quantile_method::median_unbiased)));
+        EXPECT_TRUE(xt::allclose(py_q9, xt::quantile(py_a, {0., .3, .1, 1.}, quantile_method::normal_unbiased)));
+    }
+
     /*py
     a = np.random.rand(5, 5, 5)
     */
@@ -158,59 +697,59 @@ namespace xt
     {
         // py_a
         xarray<double> py_a = {
-            {{0.0650515929852795, 0.9488855372533332, 0.9656320330745594, 0.8083973481164611, 0.3046137691733707},
-             {0.0976721140063839, 0.6842330265121569, 0.4401524937396013, 0.1220382348447788, 0.4951769101112702},
-             {0.0343885211152184, 0.9093204020787821, 0.2587799816000169, 0.662522284353982, 0.311711076089411},
-             {0.5200680211778108, 0.5467102793432796, 0.184854455525527, 0.9695846277645586, 0.7751328233611146},
-             {0.9394989415641891, 0.8948273504276488, 0.5978999788110851, 0.9218742350231168, 0.0884925020519195}},
-
-            {{0.1959828624191452, 0.0452272889105381, 0.3253303307632643, 0.388677289689482, 0.2713490317738959},
-             {0.8287375091519293, 0.3567533266935893, 0.2809345096873808, 0.5426960831582485, 0.1409242249747626},
-             {0.8021969807540397, 0.0745506436797708, 0.9868869366005173, 0.7722447692966574, 0.1987156815341724},
-             {0.0055221171236024, 0.8154614284548342, 0.7068573438476171, 0.7290071680409873, 0.7712703466859457},
-             {0.0740446517340904, 0.3584657285442726, 0.1158690595251297, 0.8631034258755935, 0.6232981268275579}},
-
-            {{0.3308980248526492, 0.0635583502860236, 0.3109823217156622, 0.325183322026747, 0.7296061783380641},
-             {0.6375574713552131, 0.8872127425763265, 0.4722149251619493, 0.1195942459383017, 0.713244787222995},
-             {0.7607850486168974, 0.5612771975694962, 0.770967179954561, 0.4937955963643907, 0.5227328293819941},
-             {0.4275410183585496, 0.0254191267440952, 0.1078914269933045, 0.0314291856867343, 0.6364104112637804},
-             {0.3143559810763267, 0.5085706911647028, 0.907566473926093, 0.2492922291488749, 0.4103829230356297}},
-
-            {{0.7555511385430487, 0.2287981654916225, 0.076979909828793, 0.289751452913768, 0.1612212872540044},
-             {0.9296976523425731, 0.808120379564417, 0.6334037565104235, 0.8714605901877177, 0.8036720768991145},
-             {0.1865700588860358, 0.8925589984899778, 0.5393422419156507, 0.8074401551640625, 0.8960912999234932},
-             {0.3180034749718639, 0.1100519245276768, 0.2279351625419417, 0.4271077886262563, 0.8180147659224931},
-             {0.8607305832563434, 0.0069521305311907, 0.5107473025775657, 0.417411003148779, 0.2221078104707302}},
-
-            {{0.1198653673336828, 0.337615171403628, 0.9429097039125192, 0.3232029320207552, 0.5187906217433661},
-             {0.7030189588951778, 0.363629602379294, 0.9717820827209607, 0.9624472949421112, 0.2517822958253642},
-             {0.4972485058923855, 0.3008783098167697, 0.2848404943774676, 0.0368869473545328, 0.6095643339798968},
-             {0.5026790232288615, 0.0514787512499894, 0.2786464642366114, 0.9082658859666537, 0.2395618906669724},
-             {0.1448948720912231, 0.489452760277563, 0.9856504541106007, 0.2420552715115004, 0.6721355474058786}}};
+            {{0.8670723185801037, 0.9132405525564713, 0.5113423988609378, 0.5015162946871996, 0.7982951789667752},
+             {0.6499639307777652, 0.7019668772577033, 0.795792669436101, 0.8900053418175663, 0.3379951568515358},
+             {0.375582952639944, 0.093981939840869, 0.578280140996174, 0.0359422737967421, 0.4655980181324602},
+             {0.5426446347075766, 0.2865412521282844, 0.5908332605690108, 0.0305002499390494, 0.0373481887492144},
+             {0.8226005606596583, 0.3601906414112629, 0.1270605126518848, 0.5222432600548044, 0.7699935530986108}},
+
+            {{0.2158210274968432, 0.6228904758190003, 0.085347464993768, 0.0516817211686077, 0.531354631568148},
+             {0.5406351216101065, 0.6374299014982066, 0.7260913337226615, 0.9758520794625346, 0.5163003483011953},
+             {0.322956472941246, 0.7951861947687037, 0.2708322512620742, 0.4389714207056361, 0.078456381342266},
+             {0.0253507434154575, 0.9626484146779251, 0.8359801205122058, 0.695974206093698, 0.4089529444142699},
+             {0.1732943200708458, 0.156437042671086, 0.2502428981645953, 0.5492266647061205, 0.7145959227000623}},
+
+            {{0.6601973767177313, 0.2799338969459428, 0.9548652806631941, 0.7378969166957685, 0.5543540525114007},
+             {0.6117207462343522, 0.4196000624277899, 0.2477309895011575, 0.3559726786512616, 0.7578461104643691},
+             {0.0143934886297559, 0.1160726405069162, 0.0460026420217527, 0.0407288023189701, 0.8554605840110072},
+             {0.7036578593800237, 0.4741738290873252, 0.0978341606510015, 0.4916158751168324, 0.4734717707805657},
+             {0.1732018699100152, 0.433851649237973, 0.3985047343973734, 0.6158500980522165, 0.6350936508676438}},
+
+            {{0.0453040097720445, 0.3746126146264712, 0.6258599157142364, 0.5031362585800877, 0.8564898411883223},
+             {0.658693631618945, 0.1629344270814297, 0.0705687474004298, 0.6424192782063156, 0.0265113105416218},
+             {0.5857755812734633, 0.9402302414249576, 0.575474177875879, 0.3881699262065219, 0.6432882184423532},
+             {0.4582528904915166, 0.5456167893159349, 0.9414648087765252, 0.3861026378007743, 0.9611905638239142},
+             {0.9053506419560637, 0.1957911347892964, 0.0693613008751655, 0.1007780013774267, 0.0182218256515497}},
+
+            {{0.0944429607559284, 0.6830067734163568, 0.071188648460229, 0.3189756302937613, 0.8448753109694546},
+             {0.0232719357358259, 0.8144684825889358, 0.2818547747733999, 0.1181648276216563, 0.6967371653641506},
+             {0.628942846779884, 0.877472013527053, 0.7350710438038858, 0.8034809303848486, 0.2820345725713065},
+             {0.1774395437797228, 0.7506147516408583, 0.806834739267264, 0.9905051420006733, 0.4126176769114265},
+             {0.3720180857927832, 0.7764129607419968, 0.3408035402530178, 0.9307573256035647, 0.8584127518430118}}};
         // py_m = np.median(a)
-        double py_m = 0.489452760277563;
+        double py_m = 0.5113423988609378;
 
         // py_m0 = np.median(a, 0)
         xarray<double> py_m0 = {
-            {0.1959828624191452, 0.2287981654916225, 0.3253303307632643, 0.325183322026747, 0.3046137691733707},
-            {0.7030189588951778, 0.6842330265121569, 0.4722149251619493, 0.5426960831582485, 0.4951769101112702},
-            {0.4972485058923855, 0.5612771975694962, 0.5393422419156507, 0.662522284353982, 0.5227328293819941},
-            {0.4275410183585496, 0.1100519245276768, 0.2279351625419417, 0.7290071680409873, 0.7712703466859457},
-            {0.3143559810763267, 0.489452760277563, 0.5978999788110851, 0.417411003148779, 0.4103829230356297}};
+            {0.2158210274968432, 0.6228904758190003, 0.5113423988609378, 0.5015162946871996, 0.7982951789667752},
+            {0.6117207462343522, 0.6374299014982066, 0.2818547747733999, 0.6424192782063156, 0.5163003483011953},
+            {0.375582952639944, 0.7951861947687037, 0.575474177875879, 0.3881699262065219, 0.4655980181324602},
+            {0.4582528904915166, 0.5456167893159349, 0.806834739267264, 0.4916158751168324, 0.4126176769114265},
+            {0.3720180857927832, 0.3601906414112629, 0.2502428981645953, 0.5492266647061205, 0.7145959227000623}};
         // py_m1 = np.median(a, 1)
         xarray<double> py_m1 = {
-            {0.0976721140063839, 0.8948273504276488, 0.4401524937396013, 0.8083973481164611, 0.311711076089411},
-            {0.1959828624191452, 0.3567533266935893, 0.3253303307632643, 0.7290071680409873, 0.2713490317738959},
-            {0.4275410183585496, 0.5085706911647028, 0.4722149251619493, 0.2492922291488749, 0.6364104112637804},
-            {0.7555511385430487, 0.2287981654916225, 0.5107473025775657, 0.4271077886262563, 0.8036720768991145},
-            {0.4972485058923855, 0.337615171403628, 0.9429097039125192, 0.3232029320207552, 0.5187906217433661}};
+            {0.6499639307777652, 0.3601906414112629, 0.578280140996174, 0.5015162946871996, 0.4655980181324602},
+            {0.2158210274968432, 0.6374299014982066, 0.2708322512620742, 0.5492266647061205, 0.5163003483011953},
+            {0.6117207462343522, 0.4196000624277899, 0.2477309895011575, 0.4916158751168324, 0.6350936508676438},
+            {0.5857755812734633, 0.3746126146264712, 0.575474177875879, 0.3881699262065219, 0.6432882184423532},
+            {0.1774395437797228, 0.7764129607419968, 0.3408035402530178, 0.8034809303848486, 0.6967371653641506}};
         // py_m2 = np.median(a, 2)
         xarray<double> py_m2 = {
-            {0.8083973481164611, 0.4401524937396013, 0.311711076089411, 0.5467102793432796, 0.8948273504276488},
-            {0.2713490317738959, 0.3567533266935893, 0.7722447692966574, 0.7290071680409873, 0.3584657285442726},
-            {0.325183322026747, 0.6375574713552131, 0.5612771975694962, 0.1078914269933045, 0.4103829230356297},
-            {0.2287981654916225, 0.808120379564417, 0.8074401551640625, 0.3180034749718639, 0.417411003148779},
-            {0.337615171403628, 0.7030189588951778, 0.3008783098167697, 0.2786464642366114, 0.489452760277563}};
+            {0.7982951789667752, 0.7019668772577033, 0.375582952639944, 0.2865412521282844, 0.5222432600548044},
+            {0.2158210274968432, 0.6374299014982066, 0.322956472941246, 0.695974206093698, 0.2502428981645953},
+            {0.6601973767177313, 0.4196000624277899, 0.0460026420217527, 0.4741738290873252, 0.433851649237973},
+            {0.5031362585800877, 0.1629344270814297, 0.5857755812734633, 0.5456167893159349, 0.1007780013774267},
+            {0.3189756302937613, 0.2818547747733999, 0.7350710438038858, 0.7506147516408583, 0.7764129607419968}};
 
         EXPECT_EQ(py_m, xt::median(py_a));
         EXPECT_EQ(py_m0, xt::median(py_a, 0));
@@ -225,125 +764,125 @@ namespace xt
     {
         // py_a
         xarray<long> py_a = {
-            {{46, 73, 109, 61, 71},
-             {33, 82, 9, 86, 90},
-             {29, 14, 114, 83, 50},
-             {4, 65, 44, 27, 12},
-             {63, 5, 57, 25, 84}},
-
-            {{21, 94, 26, 75, 43},
-             {77, 121, 124, 99, 106},
-             {123, 20, 13, 42, 3},
-             {64, 102, 40, 80, 62},
-             {78, 60, 45, 17, 28}},
-
-            {{52, 54, 24, 107, 34},
-             {67, 91, 120, 101, 7},
-             {48, 87, 93, 98, 105},
-             {53, 119, 110, 104, 76},
-             {116, 8, 31, 36, 118}},
-
-            {{11, 47, 88, 1, 59},
-             {112, 6, 41, 30, 22},
-             {49, 56, 113, 0, 81},
-             {38, 39, 32, 51, 70},
-             {95, 103, 100, 66, 89}},
-
-            {{18, 35, 108, 19, 2},
-             {92, 79, 117, 58, 72},
-             {96, 15, 111, 10, 85},
-             {69, 97, 115, 68, 23},
-             {37, 16, 122, 55, 74}}};
+            {{110, 67, 43, 114, 86},
+             {117, 31, 40, 46, 62},
+             {10, 78, 33, 103, 14},
+             {23, 101, 66, 91, 89},
+             {20, 123, 32, 50, 106}},
+
+            {{69, 108, 96, 64, 65},
+             {59, 55, 76, 19, 119},
+             {92, 2, 42, 25, 9},
+             {63, 79, 115, 30, 5},
+             {35, 3, 53, 90, 105}},
+
+            {{71, 21, 0, 44, 47},
+             {7, 102, 37, 36, 28},
+             {97, 1, 72, 26, 49},
+             {73, 81, 39, 109, 45},
+             {6, 116, 80, 100, 17}},
+
+            {{74, 34, 4, 13, 113},
+             {57, 41, 87, 38, 56},
+             {93, 121, 52, 84, 95},
+             {24, 118, 68, 15, 82},
+             {51, 94, 77, 27, 70}},
+
+            {{8, 75, 107, 60, 11},
+             {99, 48, 18, 58, 122},
+             {85, 120, 111, 83, 61},
+             {124, 16, 29, 104, 98},
+             {88, 22, 12, 54, 112}}};
 
         // py_p0 = np.partition(a, 2, 0)
         xarray<long> py_p0 = {
-            {{11, 35, 24, 1, 2}, {33, 6, 9, 30, 7}, {29, 14, 13, 0, 3}, {4, 39, 32, 27, 12}, {37, 5, 31, 17, 28}},
-
-            {{18, 47, 26, 19, 34},
-             {67, 79, 41, 58, 22},
-             {48, 15, 93, 10, 50},
-             {38, 65, 40, 51, 23},
-             {63, 8, 45, 25, 74}},
-
-            {{21, 54, 88, 61, 43},
-             {77, 82, 117, 86, 72},
-             {49, 20, 111, 42, 81},
-             {53, 97, 44, 68, 62},
-             {78, 16, 57, 36, 84}},
-
-            {{46, 94, 109, 107, 59},
-             {112, 91, 124, 101, 106},
-             {123, 56, 113, 83, 105},
-             {64, 102, 110, 80, 70},
-             {95, 103, 100, 66, 89}},
-
-            {{52, 73, 108, 75, 71},
-             {92, 121, 120, 99, 90},
-             {96, 87, 114, 98, 85},
-             {69, 119, 115, 104, 76},
-             {116, 60, 122, 55, 118}}};
+            {{8, 21, 0, 13, 11}, {7, 31, 18, 19, 28}, {10, 1, 33, 25, 9}, {23, 16, 29, 15, 5}, {6, 3, 12, 27, 17}},
+
+            {{69, 34, 4, 44, 47},
+             {57, 41, 37, 36, 56},
+             {85, 2, 42, 26, 14},
+             {24, 79, 39, 30, 45},
+             {20, 22, 32, 50, 70}},
+
+            {{71, 67, 43, 60, 65},
+             {59, 48, 40, 38, 62},
+             {92, 78, 52, 83, 49},
+             {63, 81, 66, 91, 82},
+             {35, 94, 53, 54, 105}},
+
+            {{74, 108, 96, 114, 113},
+             {117, 55, 87, 46, 119},
+             {93, 121, 72, 84, 95},
+             {73, 118, 68, 109, 89},
+             {51, 116, 77, 90, 106}},
+
+            {{110, 75, 107, 64, 86},
+             {99, 102, 76, 58, 122},
+             {97, 120, 111, 103, 61},
+             {124, 101, 115, 104, 98},
+             {88, 123, 80, 100, 112}}};
         // py_p1 = np.partition(a, 4, 1)
         xarray<long> py_p1 = {
-            {{4, 5, 44, 27, 12},
-             {29, 14, 57, 25, 50},
-             {33, 65, 9, 61, 71},
-             {46, 73, 109, 83, 84},
-             {63, 82, 114, 86, 90}},
-
-            {{64, 20, 13, 17, 3},
-             {21, 60, 26, 42, 28},
-             {77, 94, 40, 75, 43},
-             {78, 102, 45, 80, 62},
-             {123, 121, 124, 99, 106}},
-
-            {{48, 8, 24, 36, 76},
-             {52, 54, 31, 98, 34},
-             {53, 87, 93, 101, 7},
-             {67, 91, 110, 104, 105},
-             {116, 119, 120, 107, 118}},
-
-            {{38, 39, 32, 0, 70},
-             {11, 47, 88, 1, 59},
-             {49, 6, 41, 30, 22},
-             {95, 56, 100, 51, 81},
-             {112, 103, 113, 66, 89}},
-
-            {{18, 15, 108, 10, 23},
-             {37, 16, 111, 19, 2},
-             {69, 35, 115, 55, 72},
-             {92, 79, 117, 58, 74},
-             {96, 97, 122, 68, 85}}};
+            {{10, 31, 32, 91, 62},
+             {20, 67, 33, 50, 14},
+             {23, 78, 40, 46, 86},
+             {110, 101, 43, 103, 89},
+             {117, 123, 66, 114, 106}},
+
+            {{63, 2, 42, 30, 5},
+             {35, 3, 53, 25, 9},
+             {59, 55, 76, 19, 65},
+             {69, 79, 96, 64, 105},
+             {92, 108, 115, 90, 119}},
+
+            {{7, 1, 39, 36, 45},
+             {6, 21, 0, 26, 17},
+             {71, 81, 37, 44, 28},
+             {73, 102, 72, 100, 47},
+             {97, 116, 80, 109, 49}},
+
+            {{24, 41, 4, 15, 82},
+             {51, 34, 52, 13, 70},
+             {57, 94, 68, 27, 56},
+             {74, 118, 77, 38, 95},
+             {93, 121, 87, 84, 113}},
+
+            {{8, 16, 29, 58, 11},
+             {85, 22, 12, 54, 61},
+             {88, 48, 18, 60, 98},
+             {99, 75, 107, 83, 112},
+             {124, 120, 111, 104, 122}}};
         // py_p2 = np.partition(a, 3, 2)
         xarray<long> py_p2 = {
-            {{61, 46, 71, 73, 109},
-             {9, 33, 82, 86, 90},
-             {14, 29, 50, 83, 114},
-             {4, 12, 27, 44, 65},
-             {25, 57, 5, 63, 84}},
-
-            {{21, 26, 43, 75, 94},
-             {99, 77, 106, 121, 124},
-             {3, 13, 20, 42, 123},
-             {40, 62, 64, 80, 102},
-             {17, 28, 45, 60, 78}},
-
-            {{24, 34, 52, 54, 107},
-             {7, 67, 91, 101, 120},
-             {87, 48, 93, 98, 105},
-             {53, 76, 104, 110, 119},
-             {36, 31, 8, 116, 118}},
-
-            {{1, 11, 47, 59, 88},
-             {30, 22, 6, 41, 112},
-             {0, 49, 56, 81, 113},
-             {32, 38, 39, 51, 70},
-             {66, 89, 95, 100, 103}},
-
-            {{2, 18, 19, 35, 108},
-             {58, 72, 79, 92, 117},
-             {10, 85, 15, 96, 111},
-             {68, 23, 69, 97, 115},
-             {55, 37, 16, 74, 122}}};
+            {{67, 43, 86, 110, 114},
+             {46, 40, 31, 62, 117},
+             {10, 14, 33, 78, 103},
+             {23, 66, 89, 91, 101},
+             {20, 32, 50, 106, 123}},
+
+            {{64, 65, 69, 96, 108},
+             {19, 59, 55, 76, 119},
+             {25, 9, 2, 42, 92},
+             {30, 5, 63, 79, 115},
+             {3, 35, 53, 90, 105}},
+
+            {{44, 0, 21, 47, 71},
+             {7, 28, 36, 37, 102},
+             {26, 49, 1, 72, 97},
+             {39, 45, 73, 81, 109},
+             {6, 17, 80, 100, 116}},
+
+            {{13, 4, 34, 74, 113},
+             {38, 56, 41, 57, 87},
+             {84, 52, 93, 95, 121},
+             {15, 24, 68, 82, 118},
+             {27, 51, 70, 77, 94}},
+
+            {{8, 11, 60, 75, 107},
+             {58, 18, 48, 99, 122},
+             {83, 61, 85, 111, 120},
+             {16, 29, 98, 104, 124},
+             {54, 12, 22, 88, 112}}};
 
         auto p0 = xt::partition(py_a, 2, 0);
         auto p1 = xt::partition(py_a, 4, 1);
@@ -355,37 +894,37 @@ namespace xt
 
         // py_a0 = np.argpartition(a, 2, 0)
         xarray<long> py_a0 = {
-            {{3, 4, 2, 3, 4}, {0, 3, 0, 3, 2}, {0, 0, 1, 3, 1}, {0, 3, 3, 0, 0}, {4, 0, 2, 1, 1}},
+            {{4, 2, 2, 3, 4}, {2, 0, 4, 1, 2}, {0, 2, 0, 1, 1}, {0, 4, 4, 3, 1}, {2, 1, 4, 3, 2}},
 
-            {{4, 3, 1, 4, 2}, {2, 4, 3, 4, 3}, {2, 4, 2, 4, 0}, {3, 0, 1, 3, 4}, {0, 2, 1, 0, 4}},
+            {{1, 3, 3, 2, 2}, {3, 3, 2, 2, 3}, {4, 1, 1, 2, 0}, {3, 1, 2, 1, 2}, {0, 4, 0, 0, 3}},
 
-            {{1, 2, 3, 0, 1}, {1, 0, 4, 0, 4}, {3, 1, 4, 1, 3}, {2, 4, 0, 4, 1}, {1, 4, 0, 2, 0}},
+            {{2, 0, 0, 4, 1}, {1, 4, 0, 3, 0}, {1, 0, 3, 4, 2}, {1, 2, 0, 0, 3}, {1, 3, 1, 4, 1}},
 
-            {{0, 1, 0, 2, 3}, {3, 2, 1, 2, 1}, {1, 3, 3, 0, 2}, {1, 1, 2, 1, 3}, {3, 3, 3, 3, 3}},
+            {{3, 1, 1, 0, 3}, {0, 1, 3, 0, 1}, {3, 3, 2, 3, 3}, {2, 3, 3, 2, 0}, {3, 2, 3, 1, 0}},
 
-            {{2, 0, 4, 1, 0}, {4, 1, 2, 1, 0}, {4, 2, 0, 2, 4}, {4, 2, 4, 2, 2}, {2, 1, 4, 4, 2}}};
+            {{0, 4, 4, 1, 0}, {4, 2, 1, 4, 4}, {2, 4, 4, 0, 4}, {4, 0, 1, 4, 4}, {4, 0, 2, 2, 4}}};
         // py_a1 = np.argpartition(a, 4, 1)
         xarray<long> py_a1 = {
-            {{3, 4, 3, 3, 3}, {2, 2, 4, 4, 2}, {1, 3, 1, 0, 0}, {0, 0, 0, 2, 4}, {4, 1, 2, 1, 1}},
+            {{2, 1, 4, 3, 1}, {4, 0, 2, 4, 2}, {3, 2, 1, 1, 0}, {0, 3, 0, 2, 3}, {1, 4, 3, 0, 4}},
 
-            {{3, 2, 2, 4, 2}, {0, 4, 0, 2, 4}, {1, 0, 3, 0, 0}, {4, 3, 4, 3, 3}, {2, 1, 1, 1, 1}},
+            {{3, 2, 2, 3, 3}, {4, 4, 4, 2, 2}, {1, 1, 1, 1, 0}, {0, 3, 0, 0, 4}, {2, 0, 3, 4, 1}},
 
-            {{2, 4, 0, 4, 3}, {0, 0, 4, 2, 0}, {3, 2, 2, 1, 1}, {1, 1, 3, 3, 2}, {4, 3, 1, 0, 4}},
+            {{1, 2, 3, 1, 3}, {4, 0, 0, 2, 4}, {0, 3, 1, 0, 1}, {3, 1, 2, 4, 0}, {2, 4, 4, 3, 2}},
 
-            {{3, 3, 3, 2, 3}, {0, 0, 0, 0, 0}, {2, 1, 1, 1, 1}, {4, 2, 4, 3, 2}, {1, 4, 2, 4, 4}},
+            {{3, 1, 0, 3, 3}, {4, 0, 2, 0, 4}, {1, 4, 3, 4, 1}, {0, 3, 4, 1, 2}, {2, 2, 1, 2, 0}},
 
-            {{0, 2, 0, 2, 3}, {4, 4, 2, 0, 0}, {3, 0, 3, 4, 1}, {1, 1, 1, 1, 4}, {2, 3, 4, 3, 2}}};
+            {{0, 3, 3, 1, 0}, {2, 4, 4, 4, 2}, {4, 1, 1, 0, 3}, {1, 0, 0, 2, 4}, {3, 2, 2, 3, 1}}};
         // py_a2 = np.argpartition(a, 3, 2)
         xarray<long> py_a2 = {
-            {{3, 0, 4, 1, 2}, {2, 0, 1, 3, 4}, {1, 0, 4, 3, 2}, {0, 4, 3, 2, 1}, {3, 2, 1, 0, 4}},
+            {{1, 2, 4, 0, 3}, {3, 2, 1, 4, 0}, {0, 4, 2, 1, 3}, {0, 2, 4, 3, 1}, {0, 2, 3, 4, 1}},
 
-            {{0, 2, 4, 3, 1}, {3, 0, 4, 1, 2}, {4, 2, 1, 3, 0}, {2, 4, 0, 3, 1}, {3, 4, 2, 1, 0}},
+            {{3, 4, 0, 2, 1}, {3, 0, 1, 2, 4}, {3, 4, 1, 2, 0}, {3, 4, 0, 1, 2}, {1, 0, 2, 3, 4}},
 
-            {{2, 4, 0, 1, 3}, {4, 0, 1, 3, 2}, {1, 0, 2, 3, 4}, {0, 4, 3, 2, 1}, {3, 2, 1, 0, 4}},
+            {{3, 2, 1, 4, 0}, {0, 4, 3, 2, 1}, {3, 4, 1, 2, 0}, {2, 4, 0, 1, 3}, {0, 4, 2, 3, 1}},
 
-            {{3, 0, 1, 4, 2}, {3, 4, 1, 2, 0}, {3, 0, 1, 4, 2}, {2, 0, 1, 3, 4}, {3, 4, 0, 2, 1}},
+            {{3, 2, 1, 0, 4}, {3, 4, 1, 0, 2}, {3, 2, 0, 4, 1}, {3, 0, 2, 4, 1}, {3, 0, 4, 2, 1}},
 
-            {{4, 0, 3, 1, 2}, {3, 4, 1, 0, 2}, {3, 4, 1, 0, 2}, {3, 4, 0, 1, 2}, {3, 0, 1, 4, 2}}};
+            {{0, 4, 3, 1, 2}, {3, 2, 1, 0, 4}, {3, 4, 0, 2, 1}, {1, 2, 4, 3, 0}, {3, 2, 1, 0, 4}}};
 
         auto a0 = xt::argpartition(py_a, 2, 0);
         auto a1 = xt::argpartition(py_a, 4, 1);
@@ -403,125 +942,125 @@ namespace xt
     {
         // py_a
         xarray<long> py_a = {
-            {{55, 76, 33, 103, 60},
-             {2, 10, 9, 100, 67},
-             {75, 71, 40, 80, 108},
-             {78, 37, 89, 54, 74},
-             {0, 90, 11, 106, 73}},
-
-            {{20, 28, 72, 30, 63},
-             {81, 27, 97, 88, 18},
-             {39, 77, 68, 34, 29},
-             {70, 114, 3, 83, 4},
-             {101, 43, 69, 87, 119}},
-
-            {{93, 123, 122, 82, 51},
-             {48, 58, 61, 99, 79},
-             {50, 21, 109, 117, 95},
-             {91, 14, 13, 15, 7},
-             {64, 19, 44, 35, 56}},
-
-            {{6, 111, 12, 110, 102},
-             {49, 25, 41, 38, 47},
-             {42, 8, 105, 16, 1},
-             {26, 65, 22, 85, 46},
-             {107, 118, 115, 120, 57}},
-
-            {{53, 113, 24, 17, 66},
-             {32, 86, 31, 84, 62},
-             {96, 59, 121, 94, 52},
-             {124, 112, 92, 23, 36},
-             {104, 98, 5, 116, 45}}};
+            {{87, 38, 85, 104, 100},
+             {69, 50, 60, 108, 42},
+             {19, 113, 66, 122, 54},
+             {81, 31, 109, 111, 78},
+             {43, 93, 6, 105, 1}},
+
+            {{98, 110, 97, 3, 77},
+             {61, 44, 118, 8, 123},
+             {52, 96, 18, 39, 112},
+             {41, 36, 22, 119, 37},
+             {51, 121, 107, 88, 94}},
+
+            {{79, 47, 20, 120, 80},
+             {92, 33, 70, 82, 67},
+             {90, 58, 21, 84, 99},
+             {25, 10, 124, 17, 64},
+             {114, 4, 29, 55, 9}},
+
+            {{65, 24, 46, 68, 5},
+             {34, 45, 12, 28, 76},
+             {83, 32, 72, 16, 62},
+             {26, 63, 40, 106, 103},
+             {49, 59, 57, 102, 89}},
+
+            {{35, 91, 7, 27, 23},
+             {75, 101, 71, 115, 95},
+             {73, 11, 74, 56, 86},
+             {15, 2, 117, 53, 116},
+             {0, 30, 14, 13, 48}}};
 
         // py_p0 = np.partition(a, (1, 2), 0)
         xarray<long> py_p0 = {
-            {{6, 28, 12, 17, 51}, {2, 10, 9, 38, 18}, {39, 8, 40, 16, 1}, {26, 14, 3, 15, 4}, {0, 19, 5, 35, 45}},
-
-            {{20, 76, 24, 30, 60},
-             {32, 25, 31, 84, 47},
-             {42, 21, 68, 34, 29},
-             {70, 37, 13, 23, 7},
-             {64, 43, 11, 87, 56}},
-
-            {{53, 111, 33, 82, 63},
-             {48, 27, 41, 88, 62},
-             {50, 59, 105, 80, 52},
-             {78, 65, 22, 54, 36},
-             {101, 90, 44, 106, 57}},
-
-            {{55, 123, 122, 110, 102},
-             {49, 58, 61, 100, 67},
-             {75, 71, 109, 117, 108},
-             {91, 114, 89, 85, 46},
-             {107, 118, 115, 120, 119}},
-
-            {{93, 113, 72, 103, 66},
-             {81, 86, 97, 99, 79},
-             {96, 77, 121, 94, 95},
-             {124, 112, 92, 83, 74},
-             {104, 98, 69, 116, 73}}};
+            {{35, 24, 7, 3, 5}, {34, 33, 12, 8, 42}, {19, 11, 18, 16, 54}, {15, 2, 22, 17, 37}, {0, 4, 6, 13, 1}},
+
+            {{65, 38, 20, 27, 23},
+             {61, 44, 60, 28, 67},
+             {52, 32, 21, 39, 62},
+             {25, 10, 40, 53, 64},
+             {43, 30, 14, 55, 9}},
+
+            {{79, 47, 46, 68, 77},
+             {69, 45, 70, 82, 76},
+             {73, 58, 66, 56, 86},
+             {26, 31, 109, 106, 78},
+             {49, 59, 29, 88, 48}},
+
+            {{98, 110, 97, 120, 100},
+             {92, 50, 118, 108, 123},
+             {83, 96, 72, 122, 112},
+             {41, 63, 124, 111, 103},
+             {114, 93, 57, 102, 89}},
+
+            {{87, 91, 85, 104, 80},
+             {75, 101, 71, 115, 95},
+             {90, 113, 74, 84, 99},
+             {81, 36, 117, 119, 116},
+             {51, 121, 107, 105, 94}}};
         // py_p1 = np.partition(a, (1, 4), 1)
         xarray<long> py_p1 = {
-            {{0, 10, 9, 54, 60},
-             {2, 37, 11, 80, 67},
-             {55, 71, 33, 100, 73},
-             {75, 76, 40, 103, 74},
-             {78, 90, 89, 106, 108}},
-
-            {{20, 27, 3, 30, 4},
-             {39, 28, 68, 34, 18},
-             {70, 43, 69, 83, 29},
-             {81, 77, 72, 87, 63},
-             {101, 114, 97, 88, 119}},
-
-            {{48, 14, 13, 15, 7},
-             {50, 19, 44, 35, 51},
-             {64, 21, 61, 82, 56},
-             {91, 58, 109, 99, 79},
-             {93, 123, 122, 117, 95}},
-
-            {{6, 8, 12, 16, 1},
-             {26, 25, 22, 38, 46},
-             {42, 65, 41, 85, 47},
-             {49, 111, 105, 110, 57},
-             {107, 118, 115, 120, 102}},
-
-            {{32, 59, 5, 17, 36},
-             {53, 86, 24, 23, 45},
-             {96, 98, 31, 84, 52},
-             {104, 112, 92, 94, 62},
-             {124, 113, 121, 116, 66}}};
+            {{19, 31, 6, 104, 1},
+             {43, 38, 60, 105, 42},
+             {69, 50, 66, 108, 54},
+             {81, 93, 85, 111, 78},
+             {87, 113, 109, 122, 100}},
+
+            {{41, 36, 18, 3, 37},
+             {51, 44, 22, 8, 77},
+             {52, 96, 97, 39, 94},
+             {61, 110, 107, 88, 112},
+             {98, 121, 118, 119, 123}},
+
+            {{25, 4, 20, 17, 9},
+             {79, 10, 21, 55, 64},
+             {90, 33, 29, 82, 67},
+             {92, 47, 70, 84, 80},
+             {114, 58, 124, 120, 99}},
+
+            {{26, 24, 12, 16, 5},
+             {34, 32, 40, 28, 62},
+             {49, 45, 46, 68, 76},
+             {65, 59, 57, 102, 89},
+             {83, 63, 72, 106, 103}},
+
+            {{0, 2, 7, 13, 23},
+             {15, 11, 14, 27, 48},
+             {35, 30, 71, 53, 86},
+             {73, 91, 74, 56, 95},
+             {75, 101, 117, 115, 116}}};
         // py_p2 = np.partition(a, (1, 3), 2)
         xarray<long> py_p2 = {
-            {{33, 55, 60, 76, 103},
-             {2, 9, 10, 67, 100},
-             {40, 71, 75, 80, 108},
-             {37, 54, 74, 78, 89},
-             {0, 11, 73, 90, 106}},
-
-            {{20, 28, 30, 63, 72},
-             {18, 27, 81, 88, 97},
-             {29, 34, 39, 68, 77},
-             {3, 4, 70, 83, 114},
-             {43, 69, 87, 101, 119}},
-
-            {{51, 82, 93, 122, 123},
-             {48, 58, 61, 79, 99},
-             {21, 50, 95, 109, 117},
-             {7, 13, 14, 15, 91},
-             {19, 35, 44, 56, 64}},
-
-            {{6, 12, 102, 110, 111},
-             {25, 38, 41, 47, 49},
-             {1, 8, 16, 42, 105},
-             {22, 26, 46, 65, 85},
-             {57, 107, 115, 118, 120}},
-
-            {{17, 24, 53, 66, 113},
-             {31, 32, 62, 84, 86},
-             {52, 59, 94, 96, 121},
-             {23, 36, 92, 112, 124},
-             {5, 45, 98, 104, 116}}};
+            {{38, 85, 87, 100, 104},
+             {42, 50, 60, 69, 108},
+             {19, 54, 66, 113, 122},
+             {31, 78, 81, 109, 111},
+             {1, 6, 43, 93, 105}},
+
+            {{3, 77, 97, 98, 110},
+             {8, 44, 61, 118, 123},
+             {18, 39, 52, 96, 112},
+             {22, 36, 37, 41, 119},
+             {51, 88, 94, 107, 121}},
+
+            {{20, 47, 79, 80, 120},
+             {33, 67, 70, 82, 92},
+             {21, 58, 84, 90, 99},
+             {10, 17, 25, 64, 124},
+             {4, 9, 29, 55, 114}},
+
+            {{5, 24, 46, 65, 68},
+             {12, 28, 34, 45, 76},
+             {16, 32, 62, 72, 83},
+             {26, 40, 63, 103, 106},
+             {49, 57, 59, 89, 102}},
+
+            {{7, 23, 27, 35, 91},
+             {71, 75, 95, 101, 115},
+             {11, 56, 73, 74, 86},
+             {2, 15, 53, 116, 117},
+             {0, 13, 14, 30, 48}}};
 
         auto p0 = xt::partition(py_a, {1, 2}, 0);
         auto p1 = xt::partition(py_a, {1, 4}, 1);
@@ -537,37 +1076,37 @@ namespace xt
 
         // py_a0 = np.argpartition(a, (1, 2), 0)
         xarray<long> py_a0 = {
-            {{3, 1, 3, 4, 2}, {0, 0, 0, 3, 1}, {1, 3, 0, 3, 3}, {3, 2, 1, 2, 1}, {0, 2, 4, 2, 4}},
+            {{4, 3, 4, 1, 3}, {3, 2, 3, 1, 0}, {0, 4, 1, 3, 0}, {4, 4, 1, 2, 1}, {4, 2, 0, 4, 0}},
 
-            {{1, 0, 4, 1, 0}, {4, 3, 4, 4, 3}, {3, 2, 1, 1, 1}, {1, 0, 2, 4, 2}, {2, 1, 0, 1, 2}},
+            {{3, 0, 2, 4, 4}, {1, 1, 0, 3, 2}, {1, 3, 2, 1, 3}, {2, 2, 3, 4, 2}, {0, 4, 4, 2, 2}},
 
-            {{4, 3, 0, 2, 1}, {2, 1, 3, 1, 4}, {2, 4, 3, 0, 4}, {0, 3, 3, 0, 4}, {1, 0, 2, 0, 3}},
+            {{2, 2, 3, 3, 1}, {0, 3, 2, 2, 3}, {4, 2, 0, 4, 4}, {3, 0, 0, 3, 0}, {3, 3, 2, 1, 4}},
 
-            {{0, 2, 2, 3, 3}, {3, 2, 2, 0, 0}, {0, 0, 2, 2, 0}, {2, 1, 0, 3, 3}, {3, 3, 3, 3, 1}},
+            {{1, 1, 1, 2, 0}, {2, 0, 1, 0, 1}, {3, 1, 3, 0, 1}, {1, 3, 2, 0, 3}, {2, 0, 3, 3, 3}},
 
-            {{2, 4, 1, 0, 4}, {1, 4, 1, 2, 2}, {4, 1, 4, 4, 2}, {4, 4, 4, 1, 0}, {4, 4, 1, 4, 0}}};
+            {{0, 4, 0, 0, 2}, {4, 4, 4, 4, 4}, {2, 0, 4, 2, 2}, {0, 1, 4, 1, 4}, {1, 1, 1, 0, 1}}};
         // py_a1 = np.argpartition(a, (1, 4), 1)
         xarray<long> py_a1 = {
-            {{4, 1, 1, 3, 0}, {1, 3, 4, 2, 1}, {0, 2, 0, 1, 4}, {2, 0, 2, 0, 3}, {3, 4, 3, 4, 2}},
+            {{2, 3, 4, 0, 4}, {4, 0, 1, 4, 1}, {1, 1, 2, 1, 2}, {3, 4, 0, 3, 3}, {0, 2, 3, 2, 0}},
 
-            {{0, 1, 3, 0, 3}, {2, 0, 2, 2, 1}, {3, 4, 4, 3, 2}, {1, 2, 0, 4, 0}, {4, 3, 1, 1, 4}},
+            {{3, 3, 2, 0, 3}, {4, 1, 3, 1, 0}, {2, 2, 0, 2, 4}, {1, 0, 4, 4, 2}, {0, 4, 1, 3, 1}},
 
-            {{1, 3, 3, 3, 3}, {2, 4, 4, 4, 0}, {4, 2, 1, 0, 4}, {3, 1, 2, 1, 1}, {0, 0, 0, 2, 2}},
+            {{3, 4, 0, 3, 4}, {0, 3, 2, 4, 3}, {2, 1, 4, 1, 1}, {1, 0, 1, 2, 0}, {4, 2, 3, 0, 2}},
 
-            {{0, 2, 0, 2, 2}, {3, 1, 3, 1, 3}, {2, 3, 1, 3, 1}, {1, 0, 2, 0, 4}, {4, 4, 4, 4, 0}},
+            {{3, 0, 1, 2, 0}, {1, 2, 3, 1, 2}, {4, 1, 0, 0, 1}, {0, 4, 4, 4, 4}, {2, 3, 2, 3, 3}},
 
-            {{1, 2, 4, 0, 3}, {0, 1, 0, 3, 4}, {2, 4, 1, 1, 2}, {4, 3, 3, 2, 1}, {3, 0, 2, 4, 0}}};
+            {{4, 3, 0, 4, 0}, {3, 2, 4, 0, 4}, {0, 4, 1, 3, 2}, {2, 0, 2, 2, 1}, {1, 1, 3, 1, 3}}};
         // py_a2 = np.argpartition(a, (1, 3), 2)
         xarray<long> py_a2 = {
-            {{2, 0, 4, 1, 3}, {0, 2, 1, 4, 3}, {2, 1, 0, 3, 4}, {1, 3, 4, 0, 2}, {0, 2, 4, 1, 3}},
+            {{1, 2, 0, 4, 3}, {4, 1, 2, 0, 3}, {0, 4, 2, 1, 3}, {1, 4, 0, 2, 3}, {4, 2, 0, 1, 3}},
 
-            {{0, 1, 3, 4, 2}, {4, 1, 0, 3, 2}, {4, 3, 0, 2, 1}, {2, 4, 0, 3, 1}, {1, 2, 3, 0, 4}},
+            {{3, 4, 2, 0, 1}, {3, 1, 0, 2, 4}, {2, 3, 0, 1, 4}, {2, 1, 4, 0, 3}, {0, 3, 4, 2, 1}},
 
-            {{4, 3, 0, 2, 1}, {0, 1, 2, 4, 3}, {1, 0, 4, 2, 3}, {4, 2, 1, 3, 0}, {1, 3, 2, 4, 0}},
+            {{2, 1, 0, 4, 3}, {1, 4, 2, 3, 0}, {2, 1, 3, 0, 4}, {1, 3, 0, 4, 2}, {1, 4, 2, 3, 0}},
 
-            {{0, 2, 4, 3, 1}, {1, 3, 2, 4, 0}, {4, 1, 3, 0, 2}, {2, 0, 4, 1, 3}, {4, 0, 2, 1, 3}},
+            {{4, 1, 2, 0, 3}, {2, 3, 0, 1, 4}, {3, 1, 4, 2, 0}, {0, 2, 1, 4, 3}, {0, 2, 1, 4, 3}},
 
-            {{3, 2, 0, 4, 1}, {2, 0, 4, 3, 1}, {4, 1, 3, 0, 2}, {3, 4, 2, 1, 0}, {2, 4, 1, 0, 3}}};
+            {{2, 4, 3, 0, 1}, {2, 0, 4, 1, 3}, {1, 3, 0, 2, 4}, {1, 0, 3, 4, 2}, {0, 3, 2, 1, 4}}};
 
         auto a0 = xt::argpartition(py_a, {1, 2}, 0);
         auto a1 = xt::argpartition(py_a, {1, 4}, 1);
diff --git a/test/test_xsort.cpp b/test/test_xsort.cpp
index aafc7ae07..6b0e396f1 100644
--- a/test/test_xsort.cpp
+++ b/test/test_xsort.cpp
@@ -12,6 +12,7 @@
 #include "xtensor/xfixed.hpp"
 #include "xtensor/xinfo.hpp"
 #include "xtensor/xio.hpp"
+#include "xtensor/xmath.hpp"
 #include "xtensor/xrandom.hpp"
 #include "xtensor/xslice.hpp"
 #include "xtensor/xsort.hpp"
@@ -381,6 +382,40 @@ namespace xt
         EXPECT_TRUE(check_argpartition(b, r2, 2));
     }
 
+    TEST(xsort, quantile)
+    {
+        const xt::xtensor_fixed<double, xt::xshape<4, 2, 2>> data = {
+            {{3., 4.}, {2., 1.}},
+            {{-1., 1.}, {3., 2.}},
+            {{-9., -0.}, {-9., -3.}},
+            {{12., 12.}, {12., 12.}},
+        };
+
+        {
+            auto q = xt::quantile(data, {.3, .0, 1.});
+            const xt::xtensor_fixed<double, xt::xshape<3>> expected = {0.5, -9., 12.};
+            EXPECT_TRUE(xt::equal(q, expected)(0));
+        }
+
+        const xt::xtensor_fixed<double, xt::xshape<3, 4, 2>> expected1 = {
+            {{2.3, 1.9}, {0.2, 1.3}, {-9., -2.1}, {12., 12.}},
+            {{2., 1.}, {-1., 1.}, {-9., -3.}, {12., 12.}},
+            {{3., 4.}, {3., 2.}, {-9., 0.}, {12., 12.}},
+        };
+        {
+            auto q1 = xt::quantile(data, {.3, .0, 1.}, 1);
+            EXPECT_TRUE(xt::equal(q1, expected1)(0));
+        }
+        {
+            auto q1 = xt::quantile(xt::xtensor<float, 3>(data), {.3, .0, 1.}, 1);
+            EXPECT_TRUE(xt::equal(q1, expected1)(0));
+        }
+        {
+            auto q1 = xt::quantile(xt::xarray<float>(data), {.3, .0, 1.}, 1);
+            EXPECT_TRUE(xt::equal(q1, expected1)(0));
+        }
+    }
+
     TEST(xsort, median)
     {
         xt::xtensor<float, 2> a = {{3, 4, 2, 1}, {1, 1, 3, 2}, {9, 9, 9, 9}, {12, 12, 12, 12}, {5, 5, 5, 5}};