Find median of two sorted ranges

Question

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

Answer 1

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.

Answer 2

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).

Answer 3

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.

Answer 4

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;