sorting routine for forward_list

[David Kastrup <dak at gnu dot org>]

Hi,

I've dug through some old C code of mine, rewrote it to fit to
forward_list, documented it somewhat and test-drove it.  When compared
to the standard stdlibc++ installation on a current Ubuntu system,
sorting a large file:

dak@lola:~/src/mergesort$ wc /home/dak/Downloads/tex.web
  24981  125152 1030597 /home/dak/Downloads/tex.web

when compiled with -O was roughly 40% faster than using the stock
sorting routine which is suitably impressive given that it only uses
public interfaces.  Admittedly, I don't know the compiler options used
for Ubuntu's stdlibc++ but I'd hope they'd include optimization.

Of course, the actual advantage may depend a lot on the actual computer
and architecture, but the algorithm has good coherence in memory access
patterns.

It's "old school": there is not really much that could be done better
given the public interfaces of forward_list.  Sort order is stable.
Only list pointers are changed, elements are not moved or exchanged, so
iterators remain valid.

Turning it into a list-internal workhorse could give it a slight speed
boost, and that's pretty much required when using it on stuff like the
normal doubly-linked list since it does not make sense to maintain
correct backward links while sorting those: backward links can much more
efficiently be reconstituted afterwards in a single pass.

The question is who could be bothered with the task of processing this
contribution in a manner that will actually lead to its inclusion in
stdlibc++.

I have signed assignment papers for various GNU software already, so
going through with that procedure for GCC or libstdc++ would not be a
problem, given an actual interest.

-- 
David Kastrup

#include <limits>
#include <forward_list>

// Copyright (c) 2013  David Kastrup
//
// stable (sub-)list sorting routine
//
// This is a standard binary mergesort splitting (sub-)lists of size
// count into two lists of size (count/2) and ((count+1)/2), sorting
// them recursively, and using a standard list merge algorithm on the
// result.  The recursion is unwrapped.  There is a stack of
// previously sorted sublists.  The program state is not explicitly
// kept in a stack but rather deduced via binary arithmetic on the
// fly.
//
// The function returns an iterator to the last sorted element.
// It expects the length of the sublist to be sorted explicitly rather
// than in the form of an iterator.
//
// While it would be possible to code the mergesort open-ended
// (without knowledge of the total size until encountering the end of
// the list) sorting sublists of size 2^k only until it encounters an
// end iterator, the resulting subdivision of lists is non-optimal.
// In the extreme case of 2^n + 1 list elements, the maximal number of
// sorting comparisons when dividing into almost equally sized lists
// is (2^n+1)(n-1) + 3 while the maximal number of comparisons with
// sublists of size 2^k is n 2^n + 1 which is larger by 2^n - n - 1,
// almost the list size.
//
// Comparisons between two list elements will be more expensive than
// one step in walking the list is, so even if the size is not known
// in advance, calling the O(n) distance operator on the boundary
// iterators will be less expensive in the worst case than forcing the
// algorithm to run with sublists that are a power of 2.
//
// The control logic used here is actually the same as for the next
// larger power of 2 case: the length of the merged sublists is then
// adjusted to account for the missing list members.
//
// Worst case number of comparisons is
// count * ceil (log2 (count)) - (1 << ceil (log2 (count))) + 1
//
// Best case behavior (presorted data, either reverse or straight)
// takes about half the number of comparisons.  The average case for
// randomly ordered data takes hardly fewer comparisons than the worst
// case, the variance being less than 2.
//
// Sorting is stable: elements comparing equal don't change order.
//
// This implementation treats forward_list as an opaque type.  For
// internal use in predefined classes, it might make more sense to
// directly access links: in particular the use of splice_after (which
// needs to preserve list invariants) is less efficient than
// necessary.  This holds particularly for the four-argument form
// which is O(m) rather than O(1) because of the "last" iterator
// argument pointing beyond the last element with links that need to
// get changed.
//
// Sorting a doubly-linked list can be accomplished much more
// efficiently by just sorting the forward links and then
// reconstructing the backward links in a single O(n) pass afterwards.
// Preservation of type invariants would not allow such an approach to
// be used except as part of the list type implementation itself.

template <class Container, class Compare>
typename Container::const_iterator
sort_after (typename Container::const_iterator head,
	    Container &list,
	    typename Container::size_type count,
	    Compare comp)
{
  typename Container::size_type i,underpow,n,m,bitrev;
  int lev, kmax;
  typename Container::const_iterator
    headstack[std::numeric_limits<typename Container::size_type>::digits], tail, tail1, p1, p2;
  int sp = 0;

  // Recursion simulation is based on sorting lists of size 2 or more,
  // so the trivial cases 0 and 1 have to be weeded out in advance.
  switch (count)
    {
    case 0:
      return head;
    case 1:
      return ++head;
    }
  /* Set underpow to the largest power of 2 strictly smaller than count/2 */

  for (underpow = 1, kmax=0; underpow <= (count-1)/4; underpow *= 2, ++kmax)
    ;
  
  // "head" points to head of unseen chain at the start of each
  // iteration
  //
  // "i" is a running counter through list elements that have already
  // been seen.  For algorithmic reasons and to make sure it never
  // exceeds its number range, it counts quadruplets.
  // A "quadruplet" is reckoned in terms of rounding up "count" to the
  // next power of 2.  Consequently, depending on the value of
  // "count", quadruplets will be either composed of a mixture of
  // quadruplets and triplets, or of triplets and pairs.
  //
  // The way to determine which is which is by taking the bit reversal
  // of the running counter as a dithering pattern, add the absolute
  // number "count" to it and take a look at the most significant bits
  // of the result in order to determine the current sublist size.
  //
  // "bitrev" should be actually initialized to 0 like "i", but we
  // offset it already by "count".  Since that would make it prone to
  // arithmetic overflow, we subtract "underpow" again, making the
  // initial switch statement deliver cases that are too small by 1.

  lev = kmax;
  
  for (i=0, bitrev=count-underpow;;) {
    p1 = head;
    ++p1;
    p2 = p1;
    ++p2;
    if (!comp(*p2, *p1)) {	/* 50% */
      tail = p2;
    } else {
      list.splice_after (head, list, p1);
      tail = p1;
      p1 = p2;
    }
    // For any given count, only the cases 1 and 2, or the cases 2
    // and 3 will be run through.  This is the first run through.
    // After it finishes, the initial sublists will have been sorted.
    //
    // Case 1 caters for sublists of length 2 (that are already
    // sorted by above lines), case 2 caters for sublists of length
    // 3, and case 3 caters for sublists of length 4.  Sublists of 4
    // are simply done by sorting the next 2 elements and then
    // recursing into the list merge phase.

    switch (bitrev>>lev) {
      // Sorting two elements can just fall through.  This is the most
      // frequent case, so it would make some sense putting an
      // explicit case label for case 1: first to let the compiler
      // give it preferred treatment.  However, the only possible
      // values are 1, 2, and 3, anyway.
      // case 1:
      //   break;

      // Sorting three elements.  This sorts and merges sublists of
      // size 2 and 1.  In all other subdivisions, the first merged
      // sublist is never larger than the second, but this is more
      // convenient.
    case 2:
      // head points before sorted list with two members
      // p1 points to its first element
      // tail points before singleton list to be merged
      p2 = tail;
      ++p2;
      if (comp (*p2, *p1)) {
	list.splice_after (head, list, tail);
	break;
      }
      tail1 = p1++;
      if (comp (*p2, *p1)) {
	list.splice_after (tail1, list, tail);
	tail = p1;
	break;
      }
      tail = p2;
      break;

      // sublists of length four have their first sublist already
      // completed, so we shortcircuit into the loop.  Increasing
      // "lev" lets the switch exit at case 1 after the second sorted
      // pair, then breaks into the merge phase pretending that our
      // smallest chunk size is a pair rather than a quadruplet.
    case 3:
      headstack[sp++] = head;
      head = tail;
      ++lev;
      continue;
    }
    i++;
    // "head" points to second chain, "tail" to its end

    // every iteration ends with end of list in "tail", head of merged
    // list in "head"

    while (((i<<lev--)&underpow) == 0) {
      // Each iteration is one merge pass.  Successive merge passes
      // with previously merged lists are made until there is a gap in
      // the powers of two that make up the hypothetical list lengths

      // first let m be the total length of the merged list, then
      // subdivide this into floor (m/2) in n and ceil (m/2) in m.
      m = (bitrev >> lev) + 1;
      n = m/2;
      m -= n;
      tail1 = head;
      head = headstack[--sp];
      // From now on head is the predecessor of the rest of the first
      // chain which has length n, tail1 the predecessor of the rest
      // of the second chain which has length m
      p1 = head;
      ++p1;
      p2 = tail1;
      ++p2;
      do {
	if (comp (*p2, *p1))
	  {
	    do {
	      ++p2;
	    } while (--m && comp (*p2, *p1));
	    list.splice_after (head, list, tail1, p2);
	    if (!m)
	      {
		tail = tail1;
		break;
	      }
	  }
	head = p1++;
      } while (--n);
      head=headstack[sp];
    }
    if (i >= underpow)
      break;
    headstack[sp++] = head;
    head = tail;
    // i has carried into bit kmax-lev+1, so the bit reversal of i needs to
    // clear out all bits above lev, and set the lev bit
    bitrev -= underpow - (3<<lev);
    lev = kmax;
  }
  return tail;
}

// Test program follows
// gets input lines on stdin, sorts alphabetically, outputs to stdout
// Total time is written to stderr

#include <iostream>
#include <iomanip>
#include <functional>
#include <string>
#include <ctime>

using namespace std;

int
main (int ac, char *av[])
{
  forward_list<string> zip;
  string x;
  forward_list<string>::size_type n = 0;
  while (getline (cin, x))
    {
      ++n;
      zip.push_front (x);
    }
  struct timespec t1, t2;
  clock_gettime (CLOCK_PROCESS_CPUTIME_ID, &t1);
#if 1 // set to 0 to try standard library sort instead
  sort_after (zip.before_begin(), zip, n, less<string>());
#else
  zip.sort ();
#endif
  clock_gettime (CLOCK_PROCESS_CPUTIME_ID, &t2);
  cerr << (t2.tv_sec - t1.tv_sec - (t2.tv_nsec < t1.tv_nsec))
       << '.' << setw (9) << setfill ('0')
       << (t2.tv_nsec - t1.tv_nsec + 1000000000 * (t2.tv_nsec < t1.tv_nsec))
       << endl;
  for (auto p = zip.begin (); p != zip.end (); ++p)
    cout << *p << endl;
  return 0;
}