Implementing Bentley-McIlroy three-way partitioning using STL iterators?

In their talk "Quicksort is Optimal", Sedgewick and Bentley refer to a modified version of the quicksort partitioning step called Bentley-McIlroy three-way partitioning.开发者_开发知识库 This version of the partition step adapts gracefully to inputs containing equal keys by always pulling out copies of the pivot element from what remains, ensuring that when invoked on arrays containing duplicates the algorithm still performs well.

The C code for this partitioning step is reprinted here:

void threeWayPartition(Item a[], int l, int r)
{ 
  int i = l-1, j = r, p = l-1, q = r; Item v = a[r];
  if (r <= l) return;
  for (;;)
    { 
       while (a[++i] < v) ;
       while (v < a[--j]) if (j == l) break;
       if (i >= j) break;
       exch(a[i], a[j]);
       if (a[i] == v) { p++; exch(a[p], a[i]); }
       if (v == a[j]) { q--; exch(a[j], a[q]); }
    } 
  exch(a[i], a[r]); j = i-1; i = i+1;
  for (k = l; k < p; k++, j--) exch(a[k], a[j]);
  for (k = r-1; k > q; k--, i++) exch(a[i], a[k]);
}

I am interested in implementing this version of quicksort as an STL algorithm (just for my own edification, not as a replacement for the very fast std::sort). In order to do this, I would ideally accept as input to the algorithm a range of STL iterators defining the range to be sorted. Because quicksort does not require random access, I would hope that these iterators would be bidirectional iterators, since this would make the algorithm more general and would allow me to sort std::lists and other containers only supporting bidirectional access.

However, there's a slight problem with this. Notice that the very first line of the three-way partitioning algorithm contains this:

int i = l-1, p = l-1;

This ends up creating two integers that are before the range to be partitioned, which is fine because in the body of the loop they're incremented before they're used. However, if I replace these indices with bidirectional iterators, this code no longer has defined behavior because it backs an iterator up before the start of the range to sort.

My question is as follows - without substantially rewriting the core of the algorithm, is there a way to adapt this code to use STL-style iterators given that the algorithm begins by backing up an iterator at the start of a range? Right now, the only thoughts I've had would be introducing extra variables to "pretend" that we backed up the iterator on the first step, or decorating the iterators with special iterator adapters that allow you to back up before the beginning by just tracking how many logical steps you are before the start of the range. Neither of these seems very elegant... am I missing something? Is there a straightforward solution?

Thanks!

---

without substantially rewriting the core of the algorithm

That pretty much limits you to trying to hack your way around the boundary issue, so you'd need to use a custom iterator adapter, or wrap the iterator in a boost::optional or something similar so you know when its the first access.

What would be better would be to modify the algorithm to suit the tools at hand (thats exactly what the STL gets up to, using different algorithms for different iterator types).

I don't know if this is correct, but it describes the algorithm in a different way that doesn't require the iterator to go out of bounds.

---

Edit: Having said that, i've had a go at it. This code is untested as I don't know what the ouput should look like given an input - see the comments for details. It would only stand a chance of working for bidirectional/random access iterators.

#include <algorithm>
#include <iterator>

template <class Iterator>
void three_way_partition(Iterator begin, Iterator end)
{
    if (begin != end)
    {
        typename Iterator::value_type v = *(end - 1);

        // I can initialise it to begin here as its first use in the loop has
        // changed to post-increment (its pre-increment in your original
        // algorithm).
        Iterator i = begin;

        Iterator j = end - 1;

        // This should be begin - 1, but thats not valid, I set it to end
        // to act as a sentinal value, that way I know when im incrementing
        // p for the first time, and can set it to begin.
        Iterator p = end;

        Iterator q = end - 1;

        for (;;)
        {
            while (*(i++) < v);

            while (v < *(--j))
            {
                if (j == begin)
                {
                    break;
                }
            }

            if (std::distance(i, j) <= 0)
            {
                break;
            }

            if (*i == v)
            {
                if (p == end)
                {
                    p = begin;
                }
                else
                {
                    ++p;
                }

                std::iter_swap(p, i);
            }

            if (v == *j)
            {
                --q;
                std::iter_swap(j, q);
            }
        }

        std::iter_swap(i, end - 1);

        j = i - 1;
        i++;

        for (Iterator k = begin; k < p; ++k, --j)
        {
            std::iter_swap(k, j);
        }

        for (Iterator k = end - 2; k > q; --k, ++i)
        {
            std::iter_swap(i, k);
        }
    }
}

---

A major problem with the currently top ranked answer, is that the call to std::distance makes iteration quadratic in the worst case. A sequence with no unique keys will cause worse case behavior, which is especially unfortunate since this is exactly the case 3-way partitioning is designed to accelerate.

This implements Bentley-McIlroy 3-way partitioning in an optimal way for use with bi-directional iterators,

template <typename Bi1, typename Bi2>
  Bi2 swap_ranges_backward(Bi1 first1, Bi1 last1, Bi2 last2)
  {
        typedef typename std::reverse_iterator<Bi1> ri1;
        typedef typename std::reverse_iterator<Bi2> ri2;
        return std::swap_ranges(ri1(last1), ri1(first1), ri2(last2)).base();
  }

template <typename Bi, typename Cmp>
  std::pair<Bi, Bi>
  partition3(Bi first, Bi last,
    typename std::iterator_traits<Bi>::value_type pivot, Cmp comp)
  {
        Bi l_head = first;
        Bi l_tail = first;

        Bi r_head = last;
        Bi r_tail = last;

        while ( true )
         {
           // guarded to avoid overruns.
           //
           // @note this is necessary since ordered comparisons are
           // unavailable for bi-directional iterator types.
           while ( true )
              if (l_tail == r_head)
                 goto fixup_final;
              else if (comp(*l_tail, pivot))
                 ++l_tail;
              else
                 break;
           --r_head;
           while ( true )
              if (l_tail == r_head)
                 goto fixup_right;
              else if (comp(pivot, *r_head))
                 --r_head;
              else
                 break;

           std::iter_swap(l_tail, r_head);

           // compact equal to sequence front/back.
           if (!comp(*l_tail, pivot))
              std::iter_swap(l_tail, l_head++);
           if (!comp(pivot, *r_head))
              std::iter_swap(r_head, --r_tail);
           ++l_tail;
         }

fixup_right:
        // loop exited before chance to eval.
        if (!comp(pivot, *r_head))
           ++r_head;
fixup_final:
        // swap equal to partition point.
        if ((l_tail - l_head) <= (l_head - first))
           l_tail = std::swap_ranges(l_head, l_tail, first);
        else
           l_tail = swap_ranges_backward(first, l_head, l_tail);

        if ((r_tail - r_head) <= (last - r_tail))
           r_head = swap_ranges_backward(r_head, r_tail, last);
        else
           r_head = std::swap_ranges(r_tail, last, r_head);
        // equal range in values equal to pivot.
        return std::pair<Bi, Bi>(l_tail, r_head);
  }

Note: This has been tested using the Bentley verification suite. A nice side effect of the guarded advance, is that this function is safe for general purpose use (no restrictions on pivot or sequence length).

Example usage,

template<typename Bi, typename Cmp>
  void qsort_bi(Bi first, Bi last, Cmp comp)
  {
        auto nmemb = std::distance(first, last);
        if (nmemb <= 1)
           return;
        Bi pivot = first;
        std::advance(pivot, std::rand() % nmemb);

        std::pair<Bi, Bi> equal = partition3(first, last, *pivot, comp);
        qsort_bi(first, equal.first, comp);
        qsort_bi(equal.second, last, comp);
  }

template<typename Bi>
  void qsort_bi(Bi first, Bi last)
  {
        typedef typename std::iterator_traits<Bi>::value_type value_type;
        qsort_bi(first, last, std::less<value_type>());
  }

While the above sort may work, it illustrates a point another answer has already made, that is, bi-directional iterators and quicksort are a poor fit.

Without the ability to select a decent pivot in constant time, the performance hit makes quicksort an inferior choice. Additionally, bi-directional iterators are too general for an optimal sort on linked lists, because they cannot take advantage of the strengths of a list, such as constant time insertion and splicing. Finally, another more subtle (maybe arguable) problem, is that users expect sorts on linked lists to be stable.

My recommendation? The bottom-up, iterative mergesort used by the sgi STL. It's proven, stable, simple and fast (guaranteed n*log(n)). Unfortunately this algorithm doesn't seem to have a unique name, and I was unable to find a link to an implementation in isolation, so it's repeated here.

This is a very slick algorithm, it works in a way analogous to a binary counter (with non-empty lists equal to one). Counter holds lists of sizes 2^index (i.e. 1,2,4,8 ...). As each element (bit) is added, a carry may be initiated that will cascade into higher order lists (binary addition).

template <typename Tp>
  void msort_list(std::list<Tp>& in)
  {
        std::list<Tp> carry;
        std::list<Tp> counter[64];
        int fill = 1;

        while (!in.empty()) {
            carry.splice(carry.begin(), in, in.begin());
            int i = 0;
            for (; !counter[i].empty(); i++) {
                // merge upwards for stability.
                counter[i].merge(carry);
                counter[i].swap(carry);
            }
            counter[i].swap(carry);
            if (i == fill) ++fill;
        }
        for (int i = 1; i < fill; i++)
            counter[i].merge(counter[i-1]);
        in.swap(counter[fill-1]);
  }

Note: This version deviates from the original in a couple of ways. 1) We start fill at one instead of zero, this allows us to skip the size check and have the final swap work without otherwise affecting behavior. 2) The original inner loop condition adds i < fill, this check is extraneous (probably a holdover from a version where the counter array was dynamic).

---

Unfortunately, the expression "k < p" is illegal for bidirectional iterators (requires random access). It seems to me that is the real limitation you face. For example

if (i >= j) break; 

has to go away and be replaced with

if (i == j) break;

which means you'll need to add extra conditions in the "inner" loops to ensure that j (in particular) is not decremented too far. Net/net your constraint "without substantially rewriting" cannot be satisfied while making this algorithm run for bidirectional iterators.

---

Considering all the swaps that function does, wouldn't it be easier (and maybe more efficient) to just do the following,

template <typename For, typename Cmp>
  std::pair<For, For>
  partition_3way(For first, For last,
          typename std::iterator_traits<For>::value_type pivot, Cmp comp)
  {
        For lower = std::partition(first, last, std::bind2nd(comp, pivot));
        For upper = std::partition(lower, last,
                std::not1(std::bind1st(comp, pivot)));
        // return equal range for elements equal to pivot.
        return std::pair<For, For>(lower, upper);
  }