# Find median of two sorted ranges

Given two sorted ranges of elements, compute a median of all the elements in O(log(n+m)) time.

This solution contains only one call to std::lower_bound

#include <algorithm>
#include <iostream>
#include <iterator>
#include <vector>

template <typename It1, typename It2>
auto max_end(It1 begin1, It1 end1, It2 begin2, It2 end2) {
if (begin1 != end1 && begin2 != end2) {
return std::max(*std::prev(end1), *std::prev(end2));
}
if (begin1 != end1) {
return *std::prev(end1);
}
if (begin2 != end2) {
return *std::prev(end2);
}
// both ranges should never be empty
__builtin_unreachable();
}

template <typename It1, typename It2>
auto min_begin(It1 begin1, It1 end1, It2 begin2, It2 end2) {
if (begin1 != end1 && begin2 != end2) {
return std::min(*begin1, *begin2);
}
if (begin1 != end1) {
return *begin1;
}
if (begin2 != end2) {
return *begin2;
}
// both ranges should never be empty
__builtin_unreachable();
}

template<typename RandomAccessIterator1, typename RandomAccessIterator2>
auto median2(RandomAccessIterator1 begin1, RandomAccessIterator1 end1,
RandomAccessIterator2 begin2, RandomAccessIterator2 end2) {
static_assert(std::is_same<typename std::iterator_traits<RandomAccessIterator1>::value_type,
typename std::iterator_traits<RandomAccessIterator2>::value_type>::value,
"Value types of two ranges must be the same");
ssize_t n1 = std::distance(begin1, end1);
ssize_t n2 = std::distance(begin2, end2);
ssize_t mid = (n1 + n2) / 2;
auto it = std::lower_bound(begin1, end1, 42, [&](const auto & item, int){
ssize_t i1 = &item - &*begin1;
ssize_t i2 = mid - i1;
// make sure i2 is a valid index
if (i2 < 0) {return false;}
if (i2 > n2) {return true;}

return max_end(begin1, begin1+i1, begin2, begin2+i2)
>= min_begin(begin1+i1, end1, begin2+i2, end2);
});

ssize_t i1 = &*it - &*begin1;
ssize_t i2 = mid - i1;

if ((n1 + n2) % 2) {
return min_begin(begin1+i1, end1, begin2+i2, end2);
}

return (max_end(begin1, begin1+i1, begin2, begin2+i2) +
min_begin(begin1+i1, end1, begin2+i2, end2)) / 2;
}

template<typename Rng1, typename Rng2>
auto median2(Rng1 &&rng1, Rng2 &&rng2) {
using std::begin;
using std::end;
return median2(begin(rng1), end(rng1), begin(rng2), end(rng2));
}

int main() {
std::cout << median2(std::vector<float>{1, 2, 3, 4}, std::vector<float>{}) << '\n';
std::cout << median2(std::vector<float>{}, std::vector<float>{1, 2, 3, 4}) << '\n';
std::cout << median2(std::vector<float>{1, 2, 3, 4, 5}, std::vector<float>{}) << '\n';
std::cout << median2(std::vector<float>{}, std::vector<float>{1, 2, 3, 4, 5}) << '\n';
std::cout << median2(std::vector<float>{1, 2}, std::vector<float>{3}) << '\n';
std::cout << median2(std::vector<float>{1, 3}, std::vector<float>{2}) << '\n';
std::cout << median2(std::vector<float>{1, 2, 3, 4, 5}, std::vector<float>{6, 7, 8, 9, 10}) << '\n';
std::cout << median2(std::vector<float>{1, 2, 3, 4, 5}, std::vector<float>{6, 7, 8, 9}) << '\n';
std::cout << median2(std::vector<float>{1, 2, 3}, std::vector<float>{6, 7, 8, 9}) << '\n';
}


Live on Coliru

Solution to the same problem in Python on codereview

• If it’s unreachable just return the remaining value instead of the superfluous if-return. Nov 17, 2017 at 2:22

While this is not a thorough code review, here are some things you could improve or should pay attention to:

1. median2 is dangerous for containers which offer RandomAccessIterators, but aren't contiguous. In particular, ssize_t i1 = &item - &*begin1; invokes undefined behavior if item and begin1 do not belong to the same memory area, which they might if the underlying container is not contiguous (also, sometimes operator& is overwritten, which can cause a lot of problems on its own). To prepare against this, I would advise you to rename RandomAccessIterator[12] to ContiguousRandomAccessIterator[12] or something in that line of wording.

2. You are probably aware of this already, but your code is not really portable. For one, there is __builtin_unreachable, which is a gcc builtin. Also, ssize_t is a type from unistd.h, which is a posix header and might not be available on non-posix-conformant operating systems (also, you do currently not include this header, which is a bug). ssize_t is easily replaceable by std::ptrdiff_t in most cases, but if you care for that much of correctness, you could also write something like using difference_t = typename std::iterator_traits<RandomAccessIterator1>::difference_type; and use that instead. For __builtin_unreachable(), if the optimization value you get by using it is really that important for you, you should wrap it in some compiler identifying macro. If it is not important (which is probably the case), you could do a variety of things, for example replace it with a throw expression or just return a default value.

3. Rng[12] is a very confusing template parameter name, at least to me. I would expect it to stand for "random number generator", which is not fitting here. Choosing a different name would likely be beneficial here.

4. You should consider passing std::pairs instead of separate iterators. This makes clear that one begin- and one end-iterator always belong together, and helps cleaning up the parameters of your functions.

5. std::vector<float>{1, 2, 3, 4}. Please don't. This code is just ugly. The reason is that 1) you are relying on implicit conversion for all your values (why not write std::vector<float>{1.0f, 2.0f, 3.0f, 4.0f} instead?) and 2) there is no reason here to use float at all. Stick to integer types unless you do something that actually requires floating point features.

• 1. Completely agree. It there a way to assert for coontiguous property though? Yes, I should probably use addressof instead of &. Nov 16, 2017 at 17:42
• @IlyaPopov I don't think there is, sadly. I would just change the name and documentation to reflect that a contiguous iterator is required. Nov 16, 2017 at 17:44
• 2. Yes, agree completely. 3. This is the convention used by Ranges TS and range-v3 to denote ranges. I don't think it is confusing. 4. I just followed std library convention on that. 5.1) Just less typing 5.2) In the case of even number of items, I was supposed to return mean of two middle items. With floats, this is more clear. Nov 16, 2017 at 17:46
• @IlyaPopov About point 3: I see, I wasn't aware of that. About point 4: Standard library convention is not necessarily the best convention. I tend to follow the Cpp Core Guidelines, which suggest putting parameters that belong together in a pair or tuple (imo this also makes a lot of sense). About 5.1: That is not a good reason. Even if it saves you some keystrokes now, it leads to bugs with implicit conversion down the line. About 5.2: Might be, but then you shouldn't represent that by changing the element type of the vector. Either fix the output type of your function to float or... Nov 16, 2017 at 17:52
• @IlyaPopov ...have the output type as a template parameter. What you're currently doing is severely hindering usability (you can't always choose the type of your container elements freely). Nov 16, 2017 at 17:53

# Correctness

&item - &*begin1 sets you on the fast track for undefined behavior, as the value returned by *begin1 might not be inside the same allocation as item (e.g. in case of a linked-list based std::deque or an iterator returning a proxy object, as does std::vector<bool>). Operations such as comparison, addition or subtraction of pointers are only defined if they both point into the same allocation (or the first element after it) and they are both of the same type!

# Implementation

• The RandomAccessIterator property is only weakly asserted: operator+ might be implemented for a ForwardIterator (it just isn't required), so the current code would accept those in some cases, which would put the runtime above $\mathcal{O}(\log{}(n+m))$. This could be fixed with a static_assert or two.
• Is there any specific reason why rng1 and rng2 are taken as forwarding reference, as opposed to a simple const &?
• __builtin_unreachable() is a compiler extension that isn't portable.
• return (max_end(...) + min_begin(...)) / 2; might overflow (more likely for integers).
• Yes, the solution is valid only for contiguous ranges. Is there a way to assert for that? Nov 16, 2017 at 17:30
• Indeed, I should assert for RandomAccessIterators. I take ranges by forwarding ref, because this is how Ranges TS and range-v3 do it. And, you are of course correct about __builtin_unreachable, the code can easily be rewritten without it. Nov 16, 2017 at 17:40
• @IlyaPopov: No way that I know of (AFAIK they didn't introduce a contiguous_iterator_tag). That said, there are other ways to fix that (e.g. roll you own lower_bound so you have explicit knowledge about which index item is at, since that was they original intention for that "broken" piece of code). Nov 16, 2017 at 17:53

### Bug

I don't think your algorithm is correct. For example, try these two test cases:

std::cout << median2(std::vector<float>{1, 2, 3, 4, 5, 6, 7, 8}, std::vector<float>{9, 10}) << '\n';
std::cout << median2(std::vector<float>{9, 10}, std::vector<float>{1, 2, 3, 4, 5, 6, 7, 8}) << '\n';


The first results in 5.5 which is correct. The second, which is an identical case except the two vectors are switched, gives 7 which is wrong.

Other examples of incorrect test cases:

std::cout << median2(std::vector<float>{3, 4, 5, 6}, std::vector<float>{1, 2}) << '\n';
std::cout << median2(std::vector<float>{1, 4, 5, 6}, std::vector<float>{2, 3}) << '\n';


These give 3 instead of 3.5.

This test case is particularly bad:

std::cout << median2(std::vector<float>{1, 4, 5, 6}, std::vector<float>{2, 3, 7}) << '\n';


This gives me 2.66247e-44 instead of 4.

I didn't follow through all your logic but it seems to me that you are binary searching while keeping a single range across both vectors, when you should be keeping two separate ranges, one for each vector. This is based off of this solution which does it that way.

• Thank you for testing! I'll have to correct the solution. Nov 16, 2017 at 19:48
ssize_t mid = (n1 + n2) / 2;

if ((n1 + n2) % 2)


For extremely large inputs, these sums could overflow. Since this is signed arithmetic, that's Undefined Behaviour.

For the first, we can now use sd::midpoint; if stuck with C++17 or earlier, then n1 + (n2-n1) / 2 is a safe alternative, since we know they are both positive.

For the second, we could compare the numbers' parity. If exactly one of the values is odd, then the sum would be odd:

if (n1 % 2 != n2 % 2)


A similar (and more likely) overflow risk is here:

return (max_end(begin1, begin1+i1, begin2, begin2+i2) +
min_begin(begin1+i1, end1, begin2+i2, end2)) / 2;


If we don't have std::midpoint, we need to be extra careful with this one - use the addition form when the values have different sign, else the subtraction form:

auto m = max_end(begin1, begin1+i1, begin2, begin2+i2);
auto n = min_begin(begin1+i1, end1, begin2+i2, end2);
return (m<0) == (n<0)  ? m + (m - n) / 2  : (m + n) / 2;