查找第K小的數 BFPRT算法

出處 http://blog.csdn.net/adong76/article/details/10071297

BFPRT算法是解決從n個數中選擇第k大或第k小的數這個經典問題的著名算法，但很多人並不了解其細節。本文將首先介紹求解這個第k小數字問題的幾個思路，然后重點介紹在最壞情況下復雜度仍然為O(n)的BFPRT算法。

一基本思路

關於選擇第k小的數有許多方法

將n個數排序(比如快速排序或歸並排序)，選取排序后的第k個數，時間復雜度為O(nlogn)。
維護一個k個元素的最大堆，存儲當前遇到的最小的k個數，時間復雜度為O(nlogk)。這種方法同樣適用於海量數據的處理。
部分的快速排序（快速選擇算法），每次划分之后判斷第k個數在左右哪個部分，然后遞歸對應的部分，平均時間復雜度為O(n)。但最壞情況下復雜度為O(n^2)。
BFPRT算法，修改快速選擇算法的主元選取規則，使用中位數的中位數的作為主元，最壞情況下時間復雜度為O(n)。

二快速選擇算法

快速選擇算法就是修改之后的快速排序算法，前面快速排序的實現與應用這篇文章中講了它的原理和實現。

其主要思想就是在快速排序中得到划分結果之后，判斷要求的第k個數是在划分結果的左邊還是右邊，然后只處理對應的那一部分，從而達到降低復雜度的效果。

在快速排序中，平均情況下數組被划分成相等的兩部分，則時間復雜度為T(n)=2*T(n/2)+O(n)，可以解得T(n)=nlogn。
在快速選擇中，平均情況下數組也是非常相等的兩部分，但是只處理其中一部分，於是T(n)=T(n/2)+O(n)，可以解得T(n)=O(n)。

但是兩者在最壞情況下的時間復雜度均為O(n^2)，出現在每次划分之后左右總有一邊為空的情況下。為了避免這個問題，需要謹慎地選取划分的主元，一般的方法有：

固定選擇首元素或尾元素作為主元。
隨機選擇一個元素作為主元。
三數取中，選擇三個數的中位數作為主元。一般是首尾數，再加中間的一個數或者隨機的一個數。

為了方便，這里把前面的代碼也放在這里。

int partition(int a[], int l, int r) //對數組a下標從l到r的元素進行划分
{
    //隨機選取一個數作為划分的基數
    int rd = l + rand() % (r-l+1); 
    swap(a[rd], a[r]);
 
    int j = l - 1; //左邊數字最右的下標
    for (int i = l; i < r; i++)
        if (a[i] <= a[r]) 
            swap(a[++j], a[i]);
    swap(a[++j], a[r]);
    return j;
}
int NthElement(int a[], int l, int r, int id) //求數組a下標l到r中的第id個數
{
    if (l == r) return a[l];        //只有一個數
    int m = partition(a, l, r), cur = m - l + 1;
    if (id == cur) return a[m];                        //剛好是第id個數
    else if(id < cur) return NthElement(a, l, m-1, id);//第id個數在左邊
    else return(a, m+1, r, id-cur);                    //第id個數在右邊
}

三 BFPRT算法

BFPRT算法，又稱為中位數的中位數算法，由5位大牛（Blum 、 Floyd 、 Pratt 、 Rivest 、 Tarjan）提出，並以他們的名字命名。參考維基上的介紹Median of medians。

算法的思想是修改快速選擇算法的主元選取方法，提高算法在最壞情況下的時間復雜度。其主要步驟為：

首先把數組按5個數為一組進行分組，最后不足5個的忽略。對每組數進行排序（如插入排序）求取其中位數。
把上一步的所有中位數移到數組的前面，對這些中位數遞歸調用BFPRT算法求得他們的中位數。
將上一步得到的中位數作為划分的主元進行整個數組的划分。
判斷第k個數在划分結果的左邊、右邊還是恰好是划分結果本身，前兩者遞歸處理，后者直接返回答案。

首先看算法的主程序，代碼如下。小於5個數的情況直接處理返回答案。否則每5個進行求取中位數並放到數組前面，遞歸調用自身求取中位數的中位數，然后用中位數作為主元進行划分。

注意這里只利用了中位數的下標，而不關心中位數的數值，目的是方便在划分函數中使用下標直接進行交換。BFPRT算法執行完畢之后可以保證我們想要的數字是排在了它真實的位置上，所以可以直接使用中位數的下標。

int BFPRT(int a[], int l, int r, int id) //求數組a下標l到r中的第id個數
{
	if (r - l + 1 <= 5) //小於等於5個數，直接排序得到結果
	{
		insertionSort(a, l, r); return a[l + id - 1];
	}

	int t = l - 1; //當前替換到前面的中位數的下標
	for (int st = l, ed; (ed = st + 4) <= r; st += 5) //每5個進行處理
	{
		insertionSort(a, st, ed); //5個數的排序
		t++; swap(a[t], a[st+2]); //將中位數替換到數組前面，便於遞歸求取中位數的中位數
	}

	int pivotId = (l + t) >> 1; //l到t的中位數的下標，作為主元的下標
	BFPRT(a, l, t, pivotId-l+1);//不關心中位數的值，保證中位數在正確的位置
	int m = partition(a, l, r, pivotId), cur = m - l + 1;
    if (id == cur) return a[m];                   //剛好是第id個數
    else if(id < cur) return BFPRT(a, l, m-1, id);//第id個數在左邊
    else return BFPRT(a, m+1, r, id-cur);         //第id個數在右邊
}

這里的划分函數與之前稍微不同，因為指定了划分主元的下標，所以參數增加了一個，並且第一步需要交換主元的位置。代碼如下：

int partition(int a[], int l, int r, int pivotId) //對數組a下標從l到r的元素進行划分
{
    //以pivotId所在元素為划分主元
	swap(a[pivotId],a[r]);
 
    int j = l - 1; //左邊數字最右的下標
    for (int i = l; i < r; i++)
        if (a[i] <= a[r]) 
            swap(a[++j], a[i]);
    swap(a[++j], a[r]);
    return j;
}

這里簡單分析一下BFPRT算法的復雜度。

划分時以5個元素為一組求取中位數，共得到n/5個中位數，再遞歸求取中位數，復雜度為T(n/5)。

得到的中位數x作為主元進行划分，在n/5個中位數中，主元x大於其中1/2*n/5=n/10的中位數，而每個中位數在其本來的5個數的小組中又大於或等於其中的3個數，所以主元x至少大於所有數中的n/10*3=3/10*n個。同理，主元x至少小於所有數中的3/10*n個。即划分之后，任意一邊的長度至少為3/10，在最壞情況下，每次選擇都選到了7/10的那一部分，則遞歸的復雜度為T(7/10*n)。

在每5個數求中位數和划分的函數中，進行若干個次線性的掃描，其時間復雜度為c*n，其中c為常數。其總的時間復雜度滿足 T(n) <= T(n/5) + T(7/10*n) + c * n。

我們假設T(n)=x*n，其中x不一定是常數（比如x可以為n的倍數，則對應的T(n)=O(n^2)）。則有 x*n <= x*n/5 + x*7/10*n + c*n，得到 x<=10*c。於是可以知道x與n無關，T(n)<=10*c*n，為線性時間復雜度算法。而這又是最壞情況下的分析，故BFPRT可以在最壞情況下以線性時間求得n個數中的第k個數。

算法復雜度也可以用樹的方式來較准確的估計（略）