一个冷门的排序算法——平滑排序

· · 算法·理论

这个算法比耐心排序还冷门，我也费了很大的劲才学会它。它就是平滑排序（Smooth Sort）。

平滑排序算法是堆排序的变体，由 Edsger Wybe Dijkstra 于1981年提出。它是一种不稳定的比较排序算法，主要用于对大量数据进行排序。它是堆排序的变种，具有较高的效率，并且其最优时间复杂度为 O(N)，在平均、最坏情况下，它的时间复杂度也保持在 O(N \log N)，但常数因子比堆排序更大。

平滑排序算法的基本思想

平滑排序算法基于一个概念叫做“平滑堆”。它使用一个特定的堆结构，该结构使得不需要让元素下沉时几乎一定会到堆的最后面。

平滑排序会建多个按后序遍历记录的堆，每个堆的大小、遍历方法都是确定的。其中，最后一个堆的堆顶最大，倒数第二个堆的堆顶第二大，以此类推。

要弄明白这个结构，首先看齐这个数列（莱恩昂多数列）：

f(x)=\begin{cases}1&x\leqslant2\\f(x-1)+f(x-2)+1&x>2\end{cases}

比如，前几项莱恩昂多数列是：1,1,3,5,9,15,25,41,67,......

这个数列与建堆有很大的联系。如果一个平滑堆有 f(x) 个数，那么这个平滑堆的左子树就有 f(x-1) 个数，这个平滑堆的右子树就有 f(x-2) 个数，这样的话，这个平滑堆的元素量就刚好是 f(x-1)\text{（左子树的元素个数）}+f(x-2)\text{（右子树的元素个数）}+1\text{（根节点）}=f(x)，这样有重大的作用。

假如我们已经算出了莱恩昂多数列的前 i 项，那么一个大小为 f(i) ，开头为 j 的平滑堆不需要额外的空间，就可以实现 O(1) 找儿子。找右儿子可直接访问 a(j+f(i)-2)，找左儿子可直接访问 a(j+f(i-1)-1)。

那么，如果给 N 个数排序，而不是 f(i) 个数，该怎么办呢？其实对于任意一个正整数 N，一定可以用若干个莱恩昂多数的和表示，因为莱恩昂多数中包含两个1，而且每一项都是奇数。所以，平滑排序建的堆是可以计算的。这就是为什么平滑排序不使用完美二叉树，为了保证一定能建出平滑堆。

建好多个堆之后，最大的数已经在最右边，那么把最大的数当成一个独立的数，不再是最后一个堆的堆顶，此时那个堆就被切分成两个堆，现在场上的堆就多一个。

然后移动堆顶的元素，使它们重新有序。（注意，只将堆顶的元素排序，而就只有最后两个数无序，所以用插入排序，在插入时先比较这个堆的两个孩子，选出较大的，再与下一个堆的堆顶比较，若更小则继续传送元素，否则留下它，并不断将它与较大的孩子交换）排序后，次大的数和次次大的数又已经在最右边。

重复这两个操作，直到所有的数归位。

这种堆结构保证了排序操作的高效性。

平滑排序算法的美中不足

平滑排序使用了不够平衡的结构，常数大，且最坏情况常数很大，实现也不易。这就是为什么我不给出平滑排序的代码。

如果你们一定要平滑排序的代码的话，我可以提供一个非原创的代码。它的逻辑可能和我讲的不太一样。

注意，这个非原创的代码和 stdc++game.h 一样只是一个头文件，不能单独运行。

无注释版本：


#ifndef Smoothsort_Included
#define Smoothsort_Included

#include <iterator>
#include <algorithm>
#include <bitset>
template <typename RandomIterator>
void Smoothsort(RandomIterator begin, RandomIterator end);
template <typename RandomIterator, typename Comparator>
void Smoothsort(RandomIterator begin, RandomIterator end, Comparator comp);
namespace smoothsort_detail {
  const size_t kNumLeonardoNumbers = 46;
  const size_t kLeonardoNumbers[kNumLeonardoNumbers] = {
    1u, 1u, 3u, 5u, 9u, 15u, 25u, 41u, 67u, 109u, 177u, 287u, 465u, 753u, 
    1219u, 1973u, 3193u, 5167u, 8361u, 13529u, 21891u, 35421u, 57313u, 92735u,
    150049u, 242785u, 392835u, 635621u, 1028457u, 1664079u, 2692537u, 
    4356617u, 7049155u, 11405773u, 18454929u, 29860703u, 48315633u, 78176337u,
    126491971u, 204668309u, 331160281u, 535828591u, 866988873u, 1402817465u,
    2269806339u, 3672623805u
  };
  struct HeapShape {
    std::bitset<kNumLeonardoNumbers> trees;
    size_t smallestTreeSize;
  };
  template <typename RandomIterator>
  RandomIterator SecondChild(RandomIterator root) {
    return root - 1;
  }
  template <typename RandomIterator>
  RandomIterator FirstChild(RandomIterator root, size_t size) {
    return SecondChild(root) - kLeonardoNumbers[size - 2];
  }
  template <typename RandomIterator, typename Comparator>
  RandomIterator LargerChild(RandomIterator root, size_t size, Comparator comp) {
    RandomIterator first  = FirstChild(root, size);
    RandomIterator second = SecondChild(root);
    return comp(*first, *second)? second : first;
  }
  template <typename RandomIterator, typename Comparator>
  void RebalanceSingleHeap(RandomIterator root, size_t size, Comparator comp) {
    while (size > 1) {
      RandomIterator first  = FirstChild(root, size);
      RandomIterator second = SecondChild(root);
      RandomIterator largerChild;
      size_t childSize;
      if (comp(*first, *second)) { 
        largerChild = second;
        childSize = size - 2; // ... and has order k - 2.
      } else {
        largerChild = first;  // First child is larger...
        childSize = size - 1; // ... and has order k - 1.
      }
      if (!comp(*root, *largerChild))
        return;
      std::iter_swap(root, largerChild);
      root = largerChild;
      size = childSize;
    }
  }
  template <typename RandomIterator, typename Comparator>
  void LeonardoHeapRectify(RandomIterator begin, RandomIterator end,
                           HeapShape shape, Comparator comp) {
    RandomIterator itr = end - 1;
    size_t lastHeapSize;
    while (true) {
      lastHeapSize = shape.smallestTreeSize;
      if (size_t(std::distance(begin, itr)) == kLeonardoNumbers[lastHeapSize] - 1)
        break;
      RandomIterator toCompare = itr;
      if (shape.smallestTreeSize > 1) {
        RandomIterator largeChild = LargerChild(itr, shape.smallestTreeSize,
                                                comp);
        if (comp(*toCompare, *largeChild))
          toCompare = largeChild;
      }
      RandomIterator priorHeap = itr - kLeonardoNumbers[lastHeapSize];
      if (!comp(*toCompare, *priorHeap))
        break;
      std::iter_swap(itr, priorHeap);
      itr = priorHeap;
      do {
        shape.trees >>= 1;
        ++shape.smallestTreeSize;
      } while (!shape.trees[0]);
    }
    RebalanceSingleHeap(itr, lastHeapSize, comp);
  }
  template <typename RandomIterator, typename Comparator>
  void LeonardoHeapAdd(RandomIterator begin, RandomIterator end,
                       RandomIterator heapEnd,
                       HeapShape& shape, Comparator comp) {
    if (!shape.trees[0]) {
      shape.trees[0] = true;
      shape.smallestTreeSize = 1;
    }
    else if (shape.trees[1] && shape.trees[0]) {
      shape.trees >>= 2;
      shape.trees[0] = true;
      shape.smallestTreeSize += 2;
    } 
    else if (shape.smallestTreeSize == 1) {
      shape.trees <<= 1;
      shape.smallestTreeSize = 0;
      shape.trees[0] = true;
    } 
    else {
      shape.trees <<= shape.smallestTreeSize - 1;
      shape.trees[0] = true;
      shape.smallestTreeSize = 1;
    }
    bool isLast = false;
    switch (shape.smallestTreeSize) {
    case 0:
      if (end + 1 == heapEnd)
        isLast = true;
      break;
    case 1:
      if (end + 1 == heapEnd || (end + 2 == heapEnd && !shape.trees[1]))
        isLast = true;
      break;
    default:
      if (size_t(std::distance(end + 1, heapEnd)) < kLeonardoNumbers[shape.smallestTreeSize - 1] + 1)
        isLast = true;
      break;
    }
    if (!isLast)
      RebalanceSingleHeap(end, shape.smallestTreeSize, comp);
    else
      LeonardoHeapRectify(begin, end + 1, shape, comp);
  }

  template <typename RandomIterator, typename Comparator>
  void LeonardoHeapRemove(RandomIterator begin, RandomIterator end,
                          HeapShape& shape, Comparator comp) {

    if (shape.smallestTreeSize <= 1) {
      do {
        shape.trees >>= 1;
        ++shape.smallestTreeSize;
      } while (shape.trees.any() && !shape.trees[0]);
      return;
    }
    const size_t heapOrder = shape.smallestTreeSize;
    shape.trees[0] = false;
    shape.trees <<= 2;
    shape.trees[1] = shape.trees[0] = true;
    shape.smallestTreeSize -= 2;
    RandomIterator leftHeap  = FirstChild(end - 1, heapOrder);
    RandomIterator rightHeap = SecondChild(end - 1);
    HeapShape allButLast = shape;
    ++allButLast.smallestTreeSize;
    allButLast.trees >>= 1;
    LeonardoHeapRectify(begin, leftHeap + 1,  allButLast, comp);
    LeonardoHeapRectify(begin, rightHeap + 1, shape, comp);
  }
}
template <typename RandomIterator, typename Comparator>
void Smoothsort(RandomIterator begin, RandomIterator end, Comparator comp) {
  if (begin == end || begin + 1 == end) return;
  smoothsort_detail::HeapShape shape;
  shape.smallestTreeSize = 0;
  for (RandomIterator itr = begin; itr != end; ++itr)
    smoothsort_detail::LeonardoHeapAdd(begin, itr, end, shape, comp);
  for (RandomIterator itr = end; itr != begin; --itr)
    smoothsort_detail::LeonardoHeapRemove(begin, itr, shape, comp);
}
template <typename RandomIterator>
void Smoothsort(RandomIterator begin, RandomIterator end) {
  Smoothsort(begin, end,
             std::less<typename std::iterator_traits<RandomIterator>::value_type>());
}

#endif

相关链接

杨树排序算法

受平滑排序启发的一种更简单的算法是杨树排序。它以荷兰圩田中常见的一排排大小递减的树命名，对于大部分未排序的输入，它更快，但最好情况的时间复杂度退化成了 O(N \log N)。

杨树分类的重大变化在于，堆顶没有保持排序顺序。所以，每次在找最大数时，都会 O(\log N) 搜索堆顶，将最大的堆顶与最后一个数交换，这样每轮就只需调整一次堆。

作者还建议使用完美的二叉树而不是莱恩昂多树来进一步简化，但作用并不显著。

结语

有没有人发现了一个问题？

假如我们已经算出了莱恩昂多数列的前 i 项，那么一个大小为 f(i) ，开头为 j 的平滑堆不需要额外的空间，就可以实现 O(1) 找儿子。找右儿子可直接访问 a(j+f(i)-2)，找左儿子可直接访问 a(j+f(i-1)-1)。

也就是说，给定 f(i) ，好像（这个词暗藏玄机）需 O(1) 直接求出 f(i+1) 。该怎么办呢？

这个请大家自己讨论，把想法留在评论区，答案令人震惊。

感谢浏览。