Posted on April 5, 2018 in 算法

高性能排序算法

by LER0ever

马上明天就是紧张刺激的GCJ 2018的入围赛了，今天写一点关于算法的，顺便复习一下排序。
本文讨论的都是平均时间复杂度在O(n log n)以内的，涉及：

  • 传统的快排，堆排，归并
  • 稍微现代的混合(启发)/并行排序。

算法列表以及效率

排序算法最好时间平均时间最坏时间辅助空间
归并排序O(n log n)O(n log n)O(n log n)O(n)
堆排序O(n log n)O(n log n)O(n log n)O(1)
快速排序O(n log n)O(n log n)O(n^2)O(log n) ~ O(n)
自省排序O(n log n)O(n log n)O(n log n)-
PDQSortO(n)O(n log n )O(n log n)-
K路归并-O(nklog k) or O(nk^2)--
并行归并O((log n)^2)O((log n )^2)O((log n)^3)-

先来一点naïve sorting algorithms :) 就当复习着玩。

简单排序算法比较

1. 归并排序

归并排序是采用分治法的经典案例，建立在归并操作上，将两个已有序的数列合并，得到完全有序的数列。

原理

  1. 假设有两个有序序列，新开一个数组大小为两序列之和，用于存放合并后的序列。
  2. 两个指针分别初始指向两序列的起始值
  3. 比较两个指针指向的值，将较小的并入结果序列，并移动指针至序列下一个数
  4. 重复，直至某一个指针到达了该序列末端
  5. 将剩下的序列追加到结果序列

时空复杂度

最好时间最坏时间平均时间辅助空间
O(n log n)O(n log n)O(n log n)O(n)

2. 堆排序

堆排序建立在堆(完全二叉树)上，最大堆会确保父节点始终不小于任意子节点。 如果用数组做底层，那么节点i的左子树位置是2i+1，右子树为2i+2。 堆中最大元素一定在0/根节点的位置。

原理

  1. 构造最大堆
  2. 将目标数加入堆末
  3. 最大堆性质被破坏，使用沉降法调整最大堆(Reheapify)

Reheapify

  1. 同样适用分治法，给定两个已经满足最大堆性质的堆，以及一个可能不合适的根节点，执行递归操作
  2. 比较根节点和左右子节点，挑一个最大的，如果根节点已经是最大，那么整个堆就是最大堆
  3. 否则交换根节点和最大直接点，在该子树总继续调用Reheapify

时空复杂度

最好时间最坏时间平均时间辅助空间
O(n log n)O(n log n)O(n log n)O(1)

3. 快速排序

快速排序是一种改进的冒泡排序（类似的还有鸡尾酒排序, a.k.a 定向冒泡）。
每一趟排序可以将数组分为两个部分，左边均小于Pivot，右边都大于Pivot。
然后再分别对两个数组排序，直到整个数组有序。

原理

  1. 选一个Pivot元素P，将数组分为两个区间，左区间内的所有数都小于等于P，右区间都大于等于P
  2. 向左右区间分别递归，同样用执行步骤1.
  3. 重复处理，直至区间处理完毕

P的选择

快速排序最有名的特点可能就是它竟然能够退化到O(n^2)， 事实上快排性能的不稳定性就源于P选的好不好。

如果我的规则是每次选数组第一个元素，P=A[0]。那么如果出现一个倒序排列的数组， 快排的特征就相当于冒泡排序，时间复杂度退化到O(n^2)

比赛中的P选择方法
  1. 通常竞赛中的处理方法是随机化选择P，取个A[rand(0…n-1)]作为P处理。 好处就是不怎会退化，坏处就是如果数组特别大，随机数生成开销很大的情况下速度会受很大影响。
  2. 还有就是三个数取中位数的方法选择P (Median of three)， 取三个数（可以前三个，或任意三个），找出中位数把它当P
  3. 在我们学校的算法课上教授还讲了一种两个维度求中位数已得到非常接近中位数的方法， 不过我还没有实际测试性能。
曲线规避P的问题
  1. DualPivotQuickSort，取两个Pivot，实现快速三向切分的快速排序， 例如JDK1.8中的DualPivotQuicksort
  2. 元素数量较少的时候，快速排序性能比插入排序差， 所以元素数量低于某个阈值的时候采用插入排序

时空复杂度

最好时间最坏时间平均时间辅助空间
O(n log n)O(n^2)O(n log n)O(log n) ~ O(n)

值得注意的是，虽然快速排序和堆排序有同样的时间复杂度，但是一般情况下， 堆排序普遍比快排慢2-5倍。

混合/并行排序

介绍几种可能会用得到的混合排序算法, 两个基于快速排序的不稳定排序， 一个用不大到的K路外部归并排序以及一个稳定快速的并行化归并。

一. Introsort (std::sort) 内省排序

Introspective Sort 内省式排序是David R. Musser于1996年提出一种混合式排序算法。 大部分情况下特征和Median-of-3 QuickSort类似，但是当分割的结果有恶化倾向时， Introsort能够自我侦测并且切换到堆排序上，使之效率维持在O(n log n)， 并且速度优于从一开始就用堆排序。

std::sort中的Introsort是以上Introsort的改良版本，因为之前提到，在数据量比较小的情况下， 快速排序的递归代价太大，效率还不如插入排序，所以在GNU libstdc++中的 std::sort又被极限优化了：std::sort根据数据量大小以及是否退化将会自动选择采用 快速排序, 堆排序, 插入排序三者之一确保性能最大化。

原理

  1. 数据量足够大且足够随机时采用快速排序，此时效率为O(log n)
  2. 假设一旦某次分段后数据量小于一定的阈值，改用插入排序处理该段数据，此时效率可能达到O(n)
  3. 如果递归过程中发现递归层数过深或者分割行为有恶化倾向的时候，能够自动侦测并且换到堆排序，维持效率在O(n log n)

源码分析

代码来自《STL源码剖析》, GNU C++ 2.91

std::sort 函数主体
template <class RandomAccessIterator>
inline void sort(RandomAccessIterator first, RandomAccessIterator last) {
    if (first != last) {
        __introsort_loop(first, last, value_type(first), __lg(last - first) * 2);
        __final_insertion_sort(first, last);
    }
}

先判断大小有效，随后调用STL的Introsort实现，函数结束后调用插入排序

主函数
template <class RandomAccessIterator, class T, class Size>
void __introsort_loop(RandomAccessIterator first,
                      RandomAccessIterator last, T*,
                      Size depth_limit) {
    while (last - first > __stl_threshold) {
        if (depth_limit == 0) {
            partial_sort(first, last, last);
            return;
        }
        --depth_limit;
        RandomAccessIterator cut = __unguarded_partition
          (first, last, T(__median(*first, *(first + (last - first)/2),
                                   *(last - 1))));
        __introsort_loop(cut, last, value_type(first), depth_limit);
        last = cut;
    }
}

以上为函数的主要逻辑，可以看到和一般的快速排序仍然不同，看似并没有递归左子树。
而仔细观察可以发现在处理完右子树之后的last = cut使得下一次递归调用的时候处理了左子树， 所以通俗的来讲只对右边的子序列进行递归调用，左边的有循环调用处理。 这样的好处就是STL的实现能够节省将近一半的函数调用，在数据量大的时候差距就会更加明显。

同时可以看到函数中有个_median函数调用，就是之前提到了Median-of-3的快排，STL的做法是 取了第一个，最后一个，以及中间一个，三个数的中位数作为pivot，已经可以避免掉绝大多数恶化。

分割算法在STL的__unguarded_partition中可以看到

分割算法
template <class RandomAccessIterator, class T>
RandomAccessIterator __unguarded_partition(RandomAccessIterator first, 
                                           RandomAccessIterator last, 
                                           T pivot) {
    while (true) {
        while (*first < pivot) ++first;
        --last;
        while (pivot < *last) --last;
        if (!(first < last)) return first;
        iter_swap(first, last);
        ++first;
    }
}

这段就是通常快排的主体部分，用于用于给定P分割数组

递归深度：可以看到主函数源码中有一个depth_limit, 在降为0的时候会调用partial_sort， 也就是一个堆排序实现，源码如下

堆排部分排序
template <class RandomAccessIterator, class T, class Compare>
void __partial_sort(RandomAccessIterator first, RandomAccessIterator middle,
                    RandomAccessIterator last, T*, Compare comp) {
    make_heap(first, middle, comp);
    for (RandomAccessIterator i = middle; i < last; ++i)
        if (comp(*i, *first))
            __pop_heap(first, middle, i, T(*i), comp, distance_type(first));
    sort_heap(first, middle, comp);
}

template <class RandomAccessIterator, class Compare>
inline void partial_sort(RandomAccessIterator first,
                         RandomAccessIterator middle,
                         RandomAccessIterator last, Compare comp) {
    __partial_sort(first, middle, last, value_type(first), comp);
    }

另一个会改变算法的阈值为__stl_threshold，所谓的最小分段阈值，STL中的定义为16. 意思就是这种情况下采用不断递归来排序显然有些浪费，那么我们切换为插入排序，效率也能到O(n)。 所以STL终止快速排序，并且调用__final_insertion_sort处理最终的数组

最终的插入排序
template <class RandomAccessIterator>
void __final_insertion_sort(RandomAccessIterator first, 
                            RandomAccessIterator last) {
    if (last - first > __stl_threshold) {
        __insertion_sort(first, first + __stl_threshold);
        __unguarded_insertion_sort(first + __stl_threshold, last);
    }
    else
        __insertion_sort(first, last);
}

其中__unguarded_insertion_sort需要一定的条件：左区间内需要有最小值， 因为性能考虑去掉了越界检查，如果左区间内没有最小值那么一定会越界, 但是相比__insertion_sort，这个函数会显著地快(因为不需要检查越界)，所以 STL会尽可能能用这个函数的时候就不会去用__insertion_sort

时空复杂度

最好时间最坏时间平均时间辅助空间
O(n log n)O(n log n)O(n log n)O(log n)

二. Pattern-defeating Quicksort

PDQSort是目前单线程排序中最快的了(当然有更快的请指出)。 即使面对GNU libc++里的std::sort各种神优化，pdqsort仍然显著快于std::sort， 并且最快排序复杂度能够到O(n)

来自项目官方的跑分截图：

原理

pdqsort采用的Block Partitioning来自最近的一篇论文Block Quicksort
并且在分区恶化时会先尝试采用随机打乱元素的方式解决问题，如果分区过度恶化才会转为堆排序。
pdqsort在严格上升或下降或清一色数组时排序复杂度只有O(n)

实现

算法原作者的实现可以在github:orlp/pdfqsort里看到。

目前Rust社区的默认不稳定排序算法已经由原来的slice::sort换成了pdqsort， 并且Rust的实现原理上在某些情况甚至比pdqsort还要快。在这里 可以看到Rust的实现。

顺便扯一扯各类语言自带排序的实现：

  • C: quicksort
  • C++: introsort
  • D: introsort
  • Swift: introsort
  • Go: introsort
  • Crystal: introsort
  • Java: dual-pivot quicksort
  • Rust: pdqsort

时空复杂度

最好时间最坏时间平均时间辅助空间
O(n)O(n log n)O(n log n)O(log n)

三. k路归并

原理

通常我们所说的归并排序其实就是2路归并，不断将两个有序数组合并， 同样的，给你k个数组，每次挑一个最小的数放进结果数列中，就是所谓k路归并

一个简单的C语言实现，来源见注释

/*
 * main.c
 *
 *  Created on: 2011-11-11
 *      Author: liheyuan
 */
// From https://www.coder4.com/archives/2636

#include <stdio.h>
#include <stdlib.h>
#include <memory.h>

typedef int TYPE;

#define MIN_MAX 8888

void k_merge(TYPE** arrs, int* lens, int k)
{
    int* cnts = (int*) malloc(sizeof(int) * k);
    int* has_next = (int*) malloc(sizeof(int) * k);
    TYPE* datas = (int*) malloc(sizeof(int) * k);
    int i;
    int min_index = 0;
    TYPE min;

    //Read each k-way first into data
    for (i = 0; i < k; i++)
    {
        if (lens[i] >= 1)
        {
            datas[i] = arrs[i][0];
            cnts[i] = 1;
            has_next[i] = 1;
        }
        else
        {
            has_next[i] = 0;
        }
    }

    while (1)
    {

        //Select min in k
        min = MIN_MAX;
        min_index = 0;
        for (i = 0; i < k; i++)
        {
            if (has_next[i])
            {
                if (datas[i] < min)
                {
                    min = datas[i];
                    min_index = i;
                }
            }
        }

        if (min == MIN_MAX)
        {
            //No way has data
            break;
        }
        else
        {
            //print
            printf("%d\n", datas[min_index]);

            //Succ get one min
            if (cnts[min_index] < lens[min_index])
            {
                datas[min_index] = arrs[min_index][cnts[min_index]];
                cnts[min_index]++;
            }
            else
            {
                has_next[min_index] = 0;
            }
        }
    }

    free(cnts);
}

时空复杂度

因为K路归并通常用于多路外部排序，IO是算法的主导因素，竞赛可能也用不到这个算法。
通常情况下合并k个长度为n的数组复杂度为n*k^2

四. 并行归并

归并排序非常便于并行优化，因为递归分出的两部分互不影响。

ACM/ICPC由于限制编译条件，在使用C++的时候几乎无法使用并行化加速程序， 但是GCJ提供了使用Go语言的选项，所以可以使用Go语言自带的协程Goroutine实现并行。

原理

  • 并行递归，直到数据量足够小的时候用自带的IntroSort排序
  • 并行归并，一开始用自带Merge，之后并行merge四份子序列，直到最后一次合并为一个有序数列

简单易懂的Go语言优雅实现

// 并行排序
func parallelSort(b, c []int, p, q int) {
    // 数据量划分到足够小时调用普通排序
    // 在Go语言里应该是IntroSort
    if q-p <= 4096 {
        sort.Ints(b[p:q])
        return
    }

    n0 := p + (q-p)/4
    n1 := n0 + (q-p)/4
    n2 := n1 + (q-p)/4

    // 利用Goroutine并行排序
    var wg sync.WaitGroup
    wg.Add(4)
    go func() { parallelSort(b, c, p, n0); wg.Done() }()
    go func() { parallelSort(b, c, n0, n1); wg.Done() }()
    go func() { parallelSort(b, c, n1, n2); wg.Done() }()
    go func() { parallelSort(b, c, n2, q); wg.Done() }()
    wg.Wait()
    
    // Goroutine 并行归并
    wg.Add(2)
    go func() { parallelMerge(b[p:n0], b[n0:n1], c[p:n1]); wg.Done() }()
    go func() { parallelMerge(b[n1:n2], b[n2:q], c[n1:q]); wg.Done() }()
    wg.Wait()

    // 合并并行合并结果
    parallelMerge(c[p:n1], c[n1:q], b[p:q])
}


// 并行归并
func parallelMerge(m1, m2 []int, m []int) {
    
    if len(m2) > len(m1) {
        m1, m2 = m2, m1
    }

    // 数据单位小时调用普通归并算法
    if len(m1) <= 4096 {
        merge(m1, m2, m)
        return    
        
    }    
    x := len(m1) / 2
    y := sort.SearchInts(m2, m1[x])

    var wg sync.WaitGroup
    wg.Add(2)
    go func() { parallelMerge(m1[:x], m2[:y], m[:x+y]); wg.Done() }()
    go func() { parallelMerge(m1[x:], m2[y:], m[x+y:]); wg.Done() }()
    wg.Wait()
}

时间复杂度

O((log n)^2) (没有严格测试以及证明过)

ACM/GCJ