Quicksort template (for sorting corresponding arrays)

Question

I needed a function for sorting two corresponding arrays. So I wrote this code to sort N corresponding arrays!

The really "metaprogrammy" way to do it would be to use std::sort with a custom zip_iterator, but since zip_iterator's operator* returns a "proxy reference" type which is not Regular, and std::sort doesn't seem to be guaranteed to work with non-Regular types, I decided it would be much safer to build my own quicksort from scratch.

My intent is that flatmap_detail::sort_together(comparator, ctrs...) should work for any ctrs... which meet the requirements of a random-access container. (So for example I can assume that ctrs.begin() and ctrs.size() both exist.) I also assume that the value_types involved are all Regular, although perhaps move-only.

---

#include <utility> // for std::swap

namespace flatmap_detail {

template<class... Its>
void swap_together(size_t i, size_t j, Its... its)
{
    using std::swap;
    int dummy[] = {
        [&](){
            auto it = its + i;
            auto jt = its + j;
            using std::swap;
            swap(*it, *jt);
            return 0;
        }() ...
    };
    (void)dummy;
}

template<class Predicate, class Head, class... Rest>
size_t partition_together(Predicate& pred, size_t left, size_t right, Head head, const Rest... rest) {
    while (left < right) {
        while (left != right && pred(*(head + left))) ++left;
        while (left != right && !pred(*(head + (right-1)))) --right;
        if (left + 1 < right) {
            flatmap_detail::swap_together(left, right-1, head, rest...);
            ++left;
            --right;
        }
    }
    return right;
}

template<class Compare, class Head, class... Rest>
void sort_together(Compare& less, size_t left, size_t right, Head head, Rest... rest) {
    if (right - left == 2) {
        // I don't like this special case, but I don't want to pay for
        // "swap, partition an array of length 1, swap back" in this
        // common case. Can this special case be subsumed into the
        // general case somehow without losing performance?
        if (less(*(head + left), *(head + (left+1)))) {
            // nothing to do
        } else {
            flatmap_detail::swap_together(left, left+1, head, rest...);
        }
    } else if (right - left >= 3) {
        size_t pivot_idx = left + (right - left) / 2;
        // Swap the pivot element all the way to the right.
        if (pivot_idx != right - 1) {
            flatmap_detail::swap_together(pivot_idx, right-1, head, rest...);
        }
        const auto& pivot_elt = *(head + (right-1));
        auto less_than_pivot = [&](const auto& x) -> bool {
            return less(x, pivot_elt);
        };
        size_t correct_pivot_idx = flatmap_detail::partition_together(less_than_pivot, left, right-1, head, rest...);
        if (correct_pivot_idx != right-1) {
            flatmap_detail::swap_together(correct_pivot_idx, right-1, head, rest...);
        }
        flatmap_detail::sort_together(less, left, correct_pivot_idx, head, rest...);
        flatmap_detail::sort_together(less, correct_pivot_idx+1, right, head, rest...);
    }
}

template<class Compare, class Head, class... Rest>
void sort_together(Compare less, Head& head, Rest&... rest) {
    flatmap_detail::sort_together(less, 0, head.size(), head.begin(), rest.begin()...);
}

} // namespace flatmap_detail

---

And here's the test harness I used to try to find the bugs (there were plenty, at first, but I can't find any more now). I know the harness code's style could be bikeshedded to death, so I'm not that interested in comments on its style; but if you see a class of inputs that I'm failing to test adequately, I'd definitely like to know about it.

#include <algorithm>
#include <iostream>
#include <random>
#include <string>
#include <vector>

std::vector<int> iota_vector(int n) {
    std::vector<int> result;
    for (int i=0; i < n; ++i) result.emplace_back(i);
    return result;
}

std::vector<int> random_vector(int n) {
    std::vector<int> result;
    static std::mt19937 g;
    for (int i=0; i < n; ++i) result.emplace_back(g());
    return result;
}

std::vector<int> a, b;

void print_vector(const std::string& who, const std::vector<int>& v) {
    std::cout << who << ":";
    for (int elt : v) {
        std::cout << " " << elt;
    }
    std::cout << "\n";
}

void print_arrays(const std::string& when) {
    std::cout << when << "\n";
    print_vector("  a", a);
    print_vector("  b", b);
}

int main()
{
    for (int t=0; t < 1000; ++t) {
        for (int n=0; n <= 200; ++n) {
            auto oa = random_vector(n);
            auto ob = iota_vector(n);
            a = oa; b = ob;
            flatmap_detail::sort_together(std::less<int>(), a, b);

            if (!std::is_sorted(a.begin(), a.end(), std::less<int>())) {
                print_arrays("FAILURE!");
                exit(1);
            }
            for (int i=0; i < n; ++i) {
                if (a[i] != oa[b[i]]) {
                    print_arrays("FAILURE!");
                    exit(1);
                }
            }
        }
    }
}