BFPRT演算法:時間複雜度O(n)求第k小的數字(分治演算法+快排)
去年寫了一篇《分治演算法 求第 小元素 & )》的文章。介紹了一種對快排進行改進的演算法,可以在時間複雜度 的情況下,找到第 小的數字。
那時候,我還不知道這個演算法叫BFPRT演算法——現在知道了,還知道它又被稱為中位數的中位數演算法,它的最壞時間複雜度為 ,它是由Blum、Floyd、Pratt、Rivest、Tarjan提出,它的思想是修改快速選擇演算法的主元選取方法,提高演算法在最壞情況下的時間複雜度。
而且,我還發現了STL中有一個類似的函式——std::nth_element (位於標頭檔案<algorithm>
中):
#include <iostream>
#include <vector>
#include <algorithm>
#include <functional>
int main() {
std::vector<int> v{5, 6, 4, 3, 2, 6, 7, 9, 3};
std::nth_element(v.begin(), v.begin() + v.size()/2, v.end());
std::cout << "The median is " << v[v.size()/2] << '\n';
std::nth_element(v.begin(), v.begin()+1, v.end(), std::greater<int>());
std::cout << "The second largest element is " << v[1] << '\n';
}
The median is 5
The second largest element is 7
好了,言歸正傳。BFPRT演算法主要由兩部分組成(快排、基準選取函式)。基準選取函式也就是中位數的中位數演算法(Median of Medians algorithm)的實現,具體來說——就是將快排中基準選取策略進行了優化,改為每次儘可能的選擇中位數作為基準。
那麼是如何儘可能的選出中位數? 如果要找到一個精確的中位數,所消耗的時間代價將得不償失,而且對於快排演算法來說,只要基準儘可能的接近真正的中位數,就能形成近似的均分。我在上一篇文章中舉了個例子,這裡我再重複一遍:
假設,我們要找arr[18]
的近似中位數——其實,也就是找到數字8。(注意到,由於使用了分組,這裡產生的只是儘可能的中位數)
int arr[18] = { 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18 };
BFPRT演算法規定了,將5個元素作為一組(選出中位數)。然後,再選出中位數中的中位數…一直到,最終選出一個數字。
那麼,第一輪就是這樣(將選出的中位數放置在最前面):
//執行momedians之前
{ 1,2,3,4,5} {6,7,8,9,10} {11,12,13,14,15} {16,17,18 };
分別對每組sort
後選取小組中位數,再使用swap
將小組中位數放置在“最前面”對應位置(注意下文中*
號的標註,它表示這一步選出的小組中位數放置到哪兒去了)。
//開始momedians
{ 3*,2,1*,4,5};//sort->swap(3,1)
{ 3,8*,1,4,5} {6,7,2*,9,10};//sort->swap(8,2)
{ 3,8,13*,4,5} {6,7,2,9,10} {11,12,1*,14,15};//sort->swap(13,1)
{ 3,8,13,17*,5} {6,7,2,9,10} {11,12,1,14,15} {16,4*,18 };//sort->swap(17,4)
同樣,第二輪初始時就成如下樣子了,很顯然已經少於5個數字了:
//執行momedians之前
{ 3,8,13,17};
直接選出中位數8
。
可以將上述過程使用C++程式碼描述如下(對於5個數字的排序,使用Insert sort
可能會效率更高)。再次強調,下面這段程式碼唯一的作用,就是用來選出每次快排的基準:
//Median of Medians algorithm
int momedians(int* const low, int* const high) {
if (low >= high - 1) {
return *low;
}
//int grp_length = (int)std::sqrt(high - low);
int grp_length = 5, grp_idx = 0;
for (int* l = low, *r; l < high; l = r+1) {
r = (l + grp_length >= high) ? high : l + grp_length - 1;
std::sort(l, r); //可以使用下文的void isort(int* const low, int* const high)
std::swap(*(low+grp_idx++), *(l + (r - l) / 2));
}
return momedians(low, low + grp_idx);
}
寫到這裡已經將BFPRT演算法核心介紹完畢了。
如果想測試使用Insert sort
是否會帶來效率上提升的小夥伴,可以試試下面這段程式碼isort
,嘗試著替換文中的std::sort
排序函式。
//Insert sort
void isort(int* const low, int* const high) {
for (int* i = low + 1; i <= high; ++i) {
if (*i < *(i-1)) {
int border = *i, *j = i-1;
for ( ; border < *j; --j) {
*(j+1) = *j;
}
*(j+1) = border;
}
}
}
讓我們思維次再回到momedians
函式。從上述程式碼中可以看到grp_length = 5
是一個固定值。那麼,在BFPRT演算法中,為什麼是選5
個作為分組? 這個我也不是很明白,我也嘗試使用sqrt(陣列長度)
作為分組的長度。
有興趣的可以閱讀 ACdreamer 的一篇《BFPRT演算法原理》,他在文章結尾處,對為什麼使用5
作為固定分組長度進行了簡單說明。同時,還附有BFPRT演算法的最壞時間複雜度為
的證明。
好了,以下為完整的“BFPRT演算法:時間複雜度
求第k小的數字”程式碼。如上文所說,演算法主體功能是快排,只是在基準選取的時候使用了momedians
演算法——而不是直接取第一個數作為基準(嚴蔚敏版教材中的做法)。
#include <iostream>
#include <algorithm>
using namespace std;
//Median of Medians algorithm
int momedians(int* const low, int* const high) {
if (low >= high - 1) {
return *low;
}
int grp_length = 5, grp_idx = 0;
for (int* l = low, *r; l < high; l = r+1) {
r = (l + grp_length >= high) ? high : l + grp_length - 1;
std::sort(l, r);
std::swap(*(low+grp_idx++), *(l + (r - l) / 2));
}
return momedians(low, low + grp_idx);
}
//Quick sort
int qsort(int* const low, int* const high, int* const ptrk) {
int* l = low, *r = high;
if (l < r) {
int pivot = momedians(l, r);
while (l < r) {
while (l < r && *r >= pivot) {
--r;
}
*l = *r;
while (l < r && *l <= pivot) {
++l;
}
*r = *l;
}
*r = pivot;
}
//per qsort end, check left == right == ptrk?
return r == ptrk ? *ptrk :
(r > ptrk ?
qsort(low, r - 1, ptrk) :
qsort(r + 1, high, ptrk));
}
//Blum、Floyd、Pratt、Rivest、Tarjan
int bfprt(int* const low, int* const high, const int k = 1) {
if (low >= high || k < 1 || k > high - low) {
throw std::invalid_argument("low > high || k < 1");
}
return qsort(low, high-1, low + k -1);//[low, high)
}
int main() {
int arr[18] = { 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18 };
//cout << bfprt(&arr[0], &arr[0]+1, 1) << endl;
cout << bfprt(&arr[0], &arr[0] + 18, 18) << endl;
return 0;
}
以上的程式碼是針對int
陣列編寫的,下面再測試一下對於其他型別的支援情況(這裡直接粗暴的加上了一個模板)。有一點想強調的是,這裡的bfprt
函式引數是左閉右開區間[low,high)
,同時k
必須是從1
開始的正數。
//Blum、Floyd、Pratt、Rivest、Tarjan
template<typename T>
T bfprt(T* const low, T* const high, const int k = 1) {
if (low >= high || k < 1 || k > high - low) {
throw std::invalid_argument("low > high || k < 1");
}
return qsort(low, high-1, low + k -1);//[low, high)
}
在main
函式中對int
、char
、string
型別進行了測試。完整程式碼如下…
#include <iostream>
#include <algorithm>
using namespace std;
//Median of Medians algorithm
template<typename T>
T momedians(T* const low, T* const high) {
if (low >= high - 1) {
return *low;
}
int grp_length = 5, grp_idx = 0;
for (T* l = low, *r; l < high; l = r+1) {
r = (l + grp_length >= high) ? high : l + grp_length - 1;
std::sort(l, r);
std::swap(*(low+grp_idx++), *(l + (r - l) / 2));
}
return momedians(low, low + grp_idx);
}
//Quick sort
template<typename T>
T qsort(T* const low, T* const high, T* const ptrk) {
T* l = low, *r = high;
if (l < r) {
T pivot = momedians(l, r);
while (l < r) {
while (l < r && *r >= pivot) {
--r;
}
*l = *r;
while (l < r && *l <= pivot) {
++l;
}
*r = *l;
}
*r = pivot;
}
//per qsort end, check left == right == ptrk?
return r == ptrk ? *ptrk :
(r > ptrk ?
qsort(low, r - 1, ptrk) :
qsort(r + 1, high, ptrk));
}
//Blum、Floyd、Pratt、Rivest、Tarjan
template<typename T>
T bfprt(T* const low, T* const high, const int k = 1) {
if (low >= high || k < 1 || k > high - low) {
throw std::invalid_argument("low > high || k < 1");
}
return qsort(low, high-1, low + k -1);//[low, high)
}
int main() {
int arr[18] = { 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18 };
cout << bfprt(&arr[0], &arr[0] + 18, 18) << endl;
char str[7] = {'a','b','c','d','e','f','g'};
cout << bfprt(&str[0], &str[0] + 7, 4) << endl;
string s = "abcdefghijklmnopqrstuvwxyz";
cout << bfprt(&s[0], &s[0] + 26, 10) << endl;
return 0;
}
2018-12-25 北京 海淀
References:
[1] 為徑,分治演算法 求第k小元素
&