STL sort源碼剖析

STL的sort()算法，數據量大時采用Quick Sort，分段遞歸排序，一旦分段后的數據量小於某個門檻，為避免Quick Sort的遞歸調用帶來過大的額外負荷，就改用Insertion Sort。如果遞歸層次過深，還會改用Heap Sort。本文先分別介紹這個三個Sort，再整合分析STL sort算法(以上三種算法的綜合) -- Introspective Sorting(內省式排序)。

一、Insertion Sort

Insertion Sort是《算法導論》一開始就討論的算法。它的基本原理是：將初始序列的第一個元素作為一個有序序列，然后將剩下的N-1個元素按關鍵字大小依次插入序列，並一直保持有序。這個算法的復雜度為O(N^2)，最好情況下時間復雜度為O(N)。在數據量很少時，尤其還是在序列“幾近排序但尚未完成”時，有着很不錯的效果。

    // 默認以漸增方式排序
    template <class RandomAccessIterator>
    void __insertion_sort(RandomAccessIterator first,
        RandomAccessIterator last)
    {
        if (first == last) return;
                // --- insertion sort 外循環 ---
        for (RandomAccessIterator i = first + 1; i != last; ++i)
            __linear_insert(first, i, value_type(first));
        // 以上，[first,i) 形成一個子區間
    }

    template <class RandomAccessIterator, class T>
    inline void __linear_insert(RandomAccessIterator first,
        RandomAccessIterator last, T*)
    {
        T value = *last;      // 記錄尾元素
        if (value < *first){      // 尾比頭還小 (注意，頭端必為最小元素)
            copy_backward(first, last, last + 1);    // 將整個區間向右移一個位置
            *first = value;      // 令頭元素等於原先的尾元素值
        }
        else    // 尾不小於頭
            __unguarded_linear_insert(last, value);
    }

    template <class RandomAccessIterator, class T>
    void __unguarded_linear_insert(RandomAccessIterator last, T value)
    {
        RandomAccessIterator next = last;
        --next;

        // --- insertion sort 內循環 ---
        // 注意，一旦不再出現逆轉對(inversion),循環就可以結束了
        while (value < *next){    // 逆轉對(inversion)存在
            *last = *next;        // 調整
            last = next;        // 調整迭代器    
            --next;            // 左移一個位置
        }
        *last = value;            // value 的正確落腳處
    }

上述函數之所以命名為unguarded_x是因為，一般的Insertion Sort在內循環原本需要做兩次判斷，判斷是否相鄰兩元素是”逆轉對“，同時也判斷循環的行進是否超過邊界。但由於上述所示的源代碼會導致最小值必然在內循環子區間的邊緣，所以兩個判斷可合為一個判斷，所以稱為unguarded_。 省下一個判斷操作，在大數據量的情況下，影響還是可觀的 。

二、Quick Sort

Quick Sort是目前已知最快的排序法，平均復雜度為O(NlogN)，可是最壞情況下將達O(N^2)。

Quick Sort算法可以敘述如下。假設S代表將被處理的序列：

1、如果S的元素個數為0或1，結束。

2、取S中的任何一個元素，當做樞軸(pivot) v。

3、將S分割為L、R兩段，使L內的每一個元素都小於或等於v，R內的每一個元素都大於或等於v。

4、對L、R遞歸執行Quick Sort。

Median-of-Three（三點中值）

因為任何元素都可以當做樞軸(pivot)，為了避免元素輸入時不夠隨機帶來的惡化效應，最理想最穩當的方式就是取整個序列的投、尾、中央三個元素的中值(median)作為樞軸。這種做法稱為median-of-three partitioning。

    // 返回 a，b，c之居中者
    template <class T>
    inline const T& __median(const T& a, const T& b, const T& c)
    {
        if (a < b)
            if (b < c)        // a < b < c
                return b;
            else if (a < c)    // a < b, b >= c, a < c  -->     a < b <= c
                return c;
            else            // a < b, b >= c, a >= c    -->     c <= a < b
                return a;
        else if (a < c)        // c > a >= b
            return a;
        else if (b < c)        // a >= b, a >= c, b < c    -->   b < c <= a
            return c;
        else                // a >= b, a >= c, b >= c    -->      c<= b <= a
            return b;        
    }

Partitioning（分割）

分割方法有很多，以下敘述既簡單又有良好成效的做法。令first向尾移動，last向頭移動。當*first大於或等於pivot時停下來，當*last小於或等於pivot時也停下來，然后檢驗兩個迭代器是否交錯。未交錯則元素互相，然 后各自調整一個位置，再繼續相同行為。若交錯，則以此時first為軸將序列分為左右兩半，左邊值都小於或等於pivot，右邊都大於等於pivot

    template <class RandomAccessIterator, class T>
    RandomAccessIterator __unguarded_partition(
                                    RandomAccessIterator first,
                                    RandomAccessIterator last,
                                    T pivot)
    {
        while(true){
            while (*first < pivot) ++first;    // first 找到 >= pivot的元素就停
            --last;

            while (pivot < *last) --last;    // last 找到 <=pivot

            if (!(first < last)) return first;    // 交錯，結束循環    
        //    else
            iter_swap(first,last);                // 大小值交換
            ++first;                            // 調整
        }
    }

三、Heap Sort

STL中有一個partial_sort()算法。

    // paitial_sort的任務是找出middle - first個最小元素。
    template <class RandomAccessIterator>
    inline void partial_sort(RandomAccessIterator first,
                             RandomAccessIterator middle,
                             RandomAccessIterator last)
    {
        __partial_sort(first, middle, last, value_type(first));
    }
    template <class RandomAccessIterator,class T>
    inline void __partial_sort(RandomAccessIterator first,
                            RandomAccessIterator middle,
                            RandomAccessIterator last, T*)
    {
        make_heap(first, middle); // 默認是max-heap，即root是最大的
        for (RandomAccessIterator i = middle; i < last; ++i)
            if (*i < *first)
                __pop_heap(first, middle, i, T(*i), distance_type(first));
        sort_heap(first,middle);
    }

partial_sort的任務是找出middle-first個最小元素，因此，首先界定出區間[first,middle)，並利用make_heap()將它組織成一個max-heap，然后就可以講[middle,last)中的每一個元素拿來與max-heap的最大值比較(max-heap的最大值就在第一個元素)；如果小於該最大值，就互換位置並重新保持max-heap的狀態。如此一來，當我們走遍整個[middle,last)時，較大的元素都已經被抽離出[first,middle)，這時候再以sort_heap()將[first,middle)做一次排序。

由於篇幅有限，本文不再闡述堆的具體實現，建議海量Google。

四、IntroSort

不當的樞軸選擇，導致不當的分割，導致Quick Sort惡化為O(N^2)。David R. Musser於1996年提出一種混合式排序算法，Introspective Sorting。其行為在大部分情況下幾乎與 median-of-3 Quick Sort完全相同。但是當分割行為(partitioning)有惡化為二次行為傾向時，能自我偵測，轉而改用Heap Sort，使效率維持在O(NlogN)，又比一開始就使用Heap Sort來得好。大部分STL的sort內部其實就是用的IntroSort。

    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);
        }
                                                                                                        
    }
    // __lg()用來控制分割惡化的情況
    // 找出2^k <= n 的最大值,例：n=7得k=2; n=20得k=4
    template<class Size>
    inline Size __lg(Size n)
    {
        Size k;
        for (k = 0; n > 1; n >>= 1) 
            ++k;
        return k;
    }

        // 當元素個數為40時,__introsort_loop的最后一個參數
        // 即__lg(last-first)*2是5*2,意思是最多允許分割10層。

    const int  __stl_threshold = 16;

    template <class RandomAccessIterator, class T, class Size>
    void __introsort_loop(RandomAccessIterator first,
                    RandomAccessIterator last, T*, 
                    Size depth_limit)
    {
        while (last - first > __stl_threshold){        // > 16
            if (depth_limit == 0){                    // 至此，分割惡化
                partial_sort(first, last, last);    // 改用 heapsort
                return;
            }

            --depth_limit;
            // 以下是 median-of-3 partition，選擇一個夠好的樞軸並決定分割點
            // 分割點將落在迭代器cut身上
            RandomAccessIterator cut = __unguarded_partition
                (first, last, T(__median(*first,
                                         *(first + (last - first)/2),
                                        *(last - 1))));

            // 對右半段遞歸進行sort
            __introsort_loop(cut,last,value_type(first), depth_limit);

            last = cut;
            // 現在回到while循環中，准備對左半段遞歸進行sort
            // 這種寫法可讀性較差，效率也並沒有比較好
        }
    }

函數一開始就判斷序列大小，通過個數檢驗之后，再檢測分割層次，若分割層次超過指定值，就改用partial_sort()，即Heap sort。都通過了這些校驗之后，便進入與Quick Sort完全相同的程序。

當__introsort_loop()結束，[first,last)內有多個“元素個數少於或等於”16的子序列，每個序列有相當程序的排序，但尚未完全排序（因為元素個數一旦小於 __stl_threshold，就被中止了）。回到母函數，再進入__final_insertion_sort()：

    template <class RandomAccessIterator>
    void __final_insertion_sort(RandomAccessIterator first,
        RandomAccessIterator last)
    {
        if (last - first > __stl_threshold){   
            // > 16
            // 一、[first,first+16)進行插入排序
            // 二、調用__unguarded_insertion_sort，實質是直接進入插入排序內循環，
            //       *參見Insertion sort 源碼
            __insertion_sort(first,first + __stl_threshold);
            __unguarded_insertion_sort(first + __stl_threshold, last);
        }
        else
            __insertion_sort(first, last);
    }

    template <class RandomAccessIterator>
    inline void __unguarded_insertion_sort(RandomAccessIterator first,
        RandomAccessIterator last)
    {
        __unguarded_insertion_sort_aux(first, last, value_type(first));
    }

    template <class RandomAccessIterator, class T>

    void __unguarded_insertion_sort_aux(RandomAccessIterator first,
        RandomAccessIterator last,
        T*)
    {
        for (RandomAccessIterator i = first; i != last; ++i)
            __unguarded_linear_insert(i, T(*i));
    }

必須要看清楚的是，__final_insertion_sort()之前，經過__introsort_loop()的整個序列可以看成是一個個元素個數小於或等於16的子序列(注意子序列長度是不等的)，且這些子序列不但內部有相當程度排序，且更重要的是以子序列與子序列之間也是“遞增”的，意思是前一個子序列中的元素都是小於后一個子序列中的元素的，所以這個時候運用insertion_sort()，效率可是相當高的。

 

/--------------------------------------- 華麗的分割線--------------------------------------/

 

關於IntroSort的測試。

細看__introsort_loop()，是不是覺得他的快排的寫法很怪， 為什么不能直接這樣寫呢？

if (last - first > __stl_threshold){		// > 16
...
...
__introsort_loop(cut,last,value_type(first), depth_limit);
__introsort_loop(first,cut,value_type(first), depth_limit);

於是，我做了一次測試，結果如下圖：

測試結果發現，如果不像STL中那么寫，其實在數組還比較小時，還快那么一丁點，

並且即使數組變大，在一百萬條記錄時也只快0.3秒，唔。。也許STL更注重於大型數據吧，

不過像他那么寫，實在是犧牲了代碼的可讀性。

 

為什么是Insertion Sort，而不是Bubble Sort。

選擇排序(Selection sort)，插入排序(Insertion Sort)，冒泡排序(Bubble Sort)。這三個排序是初學者必須知道的三個基本排序方式，且他們速度都不快 -- O(N^2)。選擇排序就不說了，最好情況復雜度也得O(N^2)，且還是個不穩定的排序算法，直接淘汰。

可冒泡排序和插入排序相比較呢？

首先，他們都是穩定的排序算法，且最好情況下都是O(N^2)。那么我就來對他們的比較次數和移動元素次數做一次對比（最好情況下），如下：

插入排序：比較次數N-1，移動元素次數2N-1。

冒泡排序：比較次數N-1，無需移動元素。(注：我所說的冒泡排序在最基本的冒泡排序基礎上還利用了一下旗幟的方式，即尋訪完序列未發生數據交換時則表示排序已完成，無需再進行之后的比較與交換動作）

那么，這樣看來冒泡豈不是是更快，我可以把上述的__final_insertion_sort()函數改成一個__final_bubble_sort()，把每個子序列分別進行冒泡排序，豈不是更好？

事實上，具體實現時，我才發現這個想法錯了，因為寫這么一個__final_bubble_sort()，我沒有辦法確定每個子序列的大小，可我還是不甘心吶，就把bubble_sort()插在__introsort_loop()最后，這樣確實是每個子序列都用bubble_sort()又排序了一次，可是測試結果太慘了，由此可以看書Bubble Sort在“幾近排序但尚未完成”的情況下是沒多少改進作用的。

 

為什么不直接用Heap Sort

堆排序將所有的數據建成一個堆，最大的數據在堆頂，它不需要遞歸或者多維的暫存數組。算法最優最差都是O(NlogN)，不像快排，如果你人品夠差還能惡化到O(N^2)。當數據量非常大時（百萬數據），因為快排是使用遞歸設計算法的，還可能發出堆棧溢出錯誤呢。

那么為什么不直接用Heap Sort？或者說給一個最低元素閾值（__stl_threshold)時也給一個最大元素閾值（100W），即當元素數目超過這個值時，直接用Heap Sort，避免堆棧溢出呢？

對於第一個問題，我測試了一下，發現直接用Heap Sort，有時還沒有Quick Sort快呢，查閱《算法導論》發現，原來雖然Quick和Heap的時間復雜性是一樣的，但堆排序的常熟因子還是大些的，並且堆排序過程中重組堆其實也不是個省時的事。

 

 

VS2010版STL中的sort竟比我自己寫的快這么多？

首先，上文實現的這個Introsort是參照SGI STL寫的，於是，我斗膽在VS2010中拿他與std:sort比了比快慢。於是就隨機產生兩個百萬數據的vector用來測試。結果發現，VS中sort的速度竟是我的10倍以上的效率。頓時對微軟萌生敬意，可是當我仔細翻看源碼時.....
原來，microsoft的sort並沒有比sgi的sort快。只是在排序vector時，microsoft把vector的本質數據“萃取”出來了。
即，取消了vector在++時的邊界檢查語句，把vector::iterator當指針一般使用。所以才在對vector排序時會比我自己寫的introsort算法快那么多呢。

 

 

Over