# Find median of values without copying

For practice, I've been playing with calculating median values. This time, I wanted to make something that works well without copying the input values (perhaps because they are bulky, or of a type that's not copyable) and without reordering the container.

My implementation is split into two parts - find iterators to the central value(s), and combine the two values. This allows the advanced user to combine the central values in interesting ways if needed (extracting fields from composite objects, for example). But the simple median::value() interface is intended to be easy to use, and works for any type with a midpoint() function (found by ADL, or falling back to std::midpoint()).

The user can override the arithmetic type used for midpoint - that's useful when we want to get a fractional value rather than rounding for integer midpoint.

The principle is straightforward - an external partial sort on a parallel container of iterators, using std::nth_element, and std::max_element() where needed. This parallel operation means that we can be quite generous in the types we accept - we only need a forward range, where most implementations require a random-access range.

I'll present the tests first (using GoogleTest framework), as they show how the functions are used. They use containers of int, for simplicity, even though copying values rather than iterators is normally more efficient for these.

#include <gtest/gtest.h>
#include <array>
#include <forward_list>
#include <vector>

template<std::ranges::forward_range Container = std::vector<int>>
static void test_values(Container const& values, int first, int second)
{
auto const its = median::iterators(values);
EXPECT_EQ(*its.first, first);
EXPECT_EQ(*its.second, second);
}

TEST(Median, Empty)
{
EXPECT_THROW(median::iterators(std::array<bool,0>{}), std::domain_error);
}

TEST(Median, OneElement)
{
SCOPED_TRACE("");
test_values({100}, 100, 100);
}

TEST(Median, TwoElements)
{
SCOPED_TRACE("");
test_values({100, 200}, 100, 200);
SCOPED_TRACE("");
test_values({200, 100}, 100, 200);
}

TEST(Median, ThreeElements)
{
SCOPED_TRACE("");
test_values({1, 3, 2}, 2, 2);
}

TEST(Median, FourElements)
{
SCOPED_TRACE("");
test_values({8, 2, 6, 4}, 4, 6);
SCOPED_TRACE("");
test_values({4, 4, 6, 4}, 4, 4);
}

TEST(Median, FiveElements)
{
SCOPED_TRACE("");
test_values({8, 2, 6, 4, 0}, 4, 4);
}

TEST(Median, PlainArray)
{
SCOPED_TRACE("");
const int arr[] = { 2, 1, 3};
test_values(arr, 2, 2);
}

TEST(Median, CustomSort)
{
auto const values = std::array{20, 91, 92, 63, 54};
// sort by last digit
auto const compare = [](int a, int b){ return a % 10 < b % 10; };
EXPECT_EQ(median::value(values, compare), 92);
}

TEST(Median, Value)
{
auto const values = std::forward_list{0, 1, 2, 3};
auto const iters = median::iterators(values);

EXPECT_EQ(median::midval(iters), 1); // integer arithmetic
EXPECT_EQ(median::midval<double>(iters), 1.5);

// Same, but using convenience functions
EXPECT_EQ(median::value(values), 1);
EXPECT_EQ(median::value<double>(values), 1.5);

// And with reverse sort (causing integer std::midpoint() to round upwards)
EXPECT_EQ(median::value(values, std::greater<int>{}), 2);
EXPECT_EQ(median::value<double>(values, std::greater<int>{}), 1.5);
}

namespace test {
struct moveonly_int
{
int value;

moveonly_int(int i) : value{i} {}
moveonly_int(const moveonly_int&) = delete;
moveonly_int(moveonly_int&&) = default;
void operator=(const moveonly_int&) = delete;
moveonly_int& operator=(moveonly_int&&) = default;

bool operator<(const moveonly_int& other) const
{ return value < other.value; }
};

double midpoint(const moveonly_int& a, const moveonly_int& b)
{
double av = a.value;
double bv = b.value;
return std::midpoint(av, bv);
}
}

TEST(Median, MoveOnly)
{
std::array<test::moveonly_int, 4> values{0, 1, 2, 3};
EXPECT_EQ(median::value(values), 1.5);
}


Here's the implementation:

#include <algorithm>
#include <concepts>
#include <functional>
#include <iterator>
#include <numeric>
#include <ranges>
#include <utility>
#include <vector>

namespace median {

// Return a pair of iterators to the two median values
// If the input is of even length, an identical pair is returned
template<std::ranges::forward_range Range, typename Compare = std::less<>>
auto iterators(const Range& values, Compare compare = {})
-> std::pair<std::ranges::iterator_t<const Range>,
std::ranges::iterator_t<const Range>>
{
auto const begin = std::ranges::begin(values);
auto const end = std::ranges::end(values);
auto const size = std::distance(begin, end);

switch (size) {
case 0: throw std::domain_error("Attempting median of empty range");
case 1: return {begin, begin};
case 2:
auto a = begin;
auto b = a; ++b;
if (!compare(*a, *b)) { std::swap(a, b); }
return {a, b};
}

auto const it_cmp
= [compare](auto a, auto b){ return compare(*a, *b); };
std::vector<std::ranges::iterator_t<const Range>> iters;
iters.reserve(size);
for (auto it = begin;  it != end;  ++it) {
iters.push_back(it);
}
auto upper = iters.begin() + size / 2;
std::ranges::nth_element(iters, upper, it_cmp);
auto lower = size % 2 ? upper
: std::max_element(iters.begin(), upper, it_cmp);
return {*lower, *upper};
}

auto midval(auto iter_pair)
requires requires{ *iter_pair.first; *iter_pair.second; }
{
using std::midpoint;
auto const [a, b] = iter_pair;
return midpoint(*a, *b);
}

template<typename ArithType>
auto midval(auto iter_pair)
requires requires(ArithType v){ v = *iter_pair.first;
v = *iter_pair.second; }
{
using std::midpoint;
ArithType const a = *iter_pair.first;
ArithType const b = *iter_pair.second;
return midpoint(a, b);
}

template<typename ArithType>
auto value(const auto& values)
{
return midval<ArithType>(iterators(values));
}
template<typename ArithType>
auto value(const auto& values, auto compare)
{
return midval<ArithType>(iterators(values, compare));
}
auto value(const auto& values)
{
return midval(iterators(values));
}
auto value(const auto& values, auto compare)
{
return midval(iterators(values, compare));
}
}


I've compiled with plenty of warnings, and run the tests under Valgrind to eliminate any silly dangling-iterator problems.

Some specific concerns:

• Have I omitted any useful tests?
• Do I really need four overloads of median::value? I accept that the two midval() implementations are different enough that they are necessary.
• Is throwing std::domain_error an appropriate reaction to empty input?
• Is passing compare by value the correct choice? The standard algorithms do so, and I guess one can use a std::reference_wrapper to override that (if we're gathering execution statistics, perhaps).
• Have I missed any useful constraints on the template types?
• Anything else worthy of note.

I started off with a different implementation using std::multiset (passing the same test suite), which only needs to store only half as many iterators (but I suspect that set overheads, and the possibility of three O(log n) set operations per element in the second half, probably wipe out any space and time advantages, respectively). Here it is, anyway (we need to include <cassert> and <set> instead of <algorithm> and <vector>):

        std::multiset<std::ranges::iterator_t<const Range>, decltype(it_cmp)>
sorted(it_cmp);
auto const halfway = begin + size/2 + 1u;
for (auto it = begin;  it != halfway;  ++it) {
sorted.insert(it);
}
for (auto it = halfway; it != end;  ++it) {
auto last = sorted.end(); --last;
if (it_cmp(it, *sorted.begin())) {
// before first
sorted.erase(last);
} else if (it_cmp(it, *last)) {
// before last
sorted.erase(last);
sorted.erase(sorted.begin());
sorted.insert(it);
} else {
// after last
sorted.erase(sorted.begin());
}
}
// sorted contains one element for odd-length input, or two
// elements for even-length input.
assert(sorted.size() == 1u + !(size % 2));
auto m = sorted.begin();
auto n = sorted.end();
return {*m, *--n};


Although I'm not using this version, please do feel free to review it - I'm always looking to learn!

• Have I omitted any useful tests?

Yes. While you did test the edge case of an empty range, there are other edge cases and extreme values that you should test. For example, the median of {INT_MIN, INT_MAX} for example. Also, consider the median of something with plus and/or minus infinity, and the median of some array with a few NaNs randomly scattered in.

• Do I really need four overloads of median::value? I accept that the two midval() implementations are different enough that they are necessary.

You can easily get rid of one by having the first templated one use a default value for the comparator:

template<typename ArithType, typename Compare = std::less<>>
auto value(const auto& values, Compare compare = {})
{
return midval<ArithType>(iterators(values, compare));
}


If you could provide a default value for ArithType, you wouldn't need the untemplated overloads. However, you can't, unless you use a hack like:

template<typename ArithType = void, typename Compare = std::less<>>
auto value(const auto& values, Compare compare = std::less<>{})
{
if constexpr (std::is_same_v<ArithType, void>)
return midval(iterators(values, compare));
else
return midval<ArithType>(iterators(values, compare));
}

• Is throwing std::domain_error an appropriate reaction to empty input?

I would allow it, although one might also argue that std::invalid_argument is more appropriate.

• Is passing compare by value the correct choice? The standard algorithms do so, and I guess one can use a std::reference_wrapper to override that (if we're gathering execution statistics, perhaps).

I would copy the semantics of the STL.

• Have I missed any useful constraints on the template types?

Yes. I would require that the value type of the range and the comparator satisfy the std::strict_weak_order concept. Basically, that ensures that you can sort the values.

You could also define a concept that checks that midpoint() can be applied, so that trying to get the median of an array of std::strings will result in a better error message, but on the other hand some people might argue that that would be too specific a concept. Ideally, std::is_arithmetic would be a good choice, except it only allows built-in types.

• Anything else worthy of note.

See below.

# It might be more expensive to sort iterators than to copy the range

A std::vector of iterators might be a lot more expensive to sort than to just copy the values themselves into a new std::vector, especially if you just want the median of ints, floats or doubles. The iterator type of some containers can be quite large, and you pay for the extra indirection. On the other hand, your solution works well even if the input is a container of uncopyable and unmovable types, or if the value type itself is very large (for example, bignums).

# Consider adding a projection parameter

Using a custom midpoint() function works if ADL works. But consider that I want to get the median of std::complex<double> numbers, sorted on their absolute value. The only way to make this work with your code is to inject an overload for std::midpoint() into namespace std, which is not so nice. Having a projection parameter avoids this issue.

• Excellent - thank you! I considered switching to copying the range for small, simple types, but then ended up with a different interface (without the pair-of-iterators return type) and decided that was probably best left for a separate function, for the caller to choose. Good point also about the cases where ADL is no help. Commented Feb 12, 2022 at 21:46
• I guess there's only two reasonable ways to handle NaNs, since they are not part of a total ordering: prohibit them entirely (throw an exception) or exclude them from the calculation. The latter approach would require passing in yet another predicate, which is going to explode the interface yet more - I might have a play with that to teach myself some techniques to deal with it. Oh, there's a third option - just declare the behaviour to be undefined if there's not a total ordering. And that's the one I chose. ;-) Infinities and the four combinations of max and min I'll certainly add. Commented Feb 13, 2022 at 10:46
• Well, it grew arms and legs... Commented Feb 19, 2022 at 16:17