BFPRT算法是解决从n个数中选择第k大或第k小的数这个经典问题的著名算法，但很多人并不了解其细节。本文将首先介绍求解这个第k小数字问题的几个思路，然后重点介绍在最坏情况下复杂度仍然为O(n)的BFPRT算法。

一 基本思路

关于选择第k小的数字有许多方法，其效率和复杂度各不一样，可以根据实际情况进行选择。

1. 将n个数排序(比如快速排序或归并排序)，选取排序后的第k个数，时间复杂度为O(nlogn)。使用STL函数sort可以大大减少编码量。
2. 将方法1中的排序方法改为线性时间排序算法(如基数排序或计数排序)，时间复杂度为O(n)。但线性时间排序算法使用限制较多，不常使用。
3. 维护一个k个元素的最大堆，存储当前遇到的最小的k个数，时间复杂度为O(nlogk)。这种方法同样适用于海量数据的处理。
4. 部分的选择排序，即把最小的放在第1位，第二小的放在第2位，直到第k位为止，时间复杂度为O(kn)。实现非常简单。
5. 部分的快速排序（快速选择算法），每次划分之后判断第k个数在左右哪个部分，然后递归对应的部分，平均时间复杂度为O(n)。但最坏情况下复杂度为O(n^2)。
6. 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)，出现在每次划分之后左右总有一边为空的情况下。为了避免这个问题，需要谨慎地选取划分的主元，一般的方法有：

1. 固定选择首元素或尾元素作为主元。
2. 随机选择一个元素作为主元。
3. 三数取中，选择三个数的中位数作为主元。一般是首尾数，再加中间的一个数或者随机的一个数。

为了方便，这里把前面的代码也放在这里。

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。

算法的思想是修改快速选择算法的主元选取方法，提高算法在最坏情况下的时间复杂度。其主要步骤为：

1. 首先把数组按5个数为一组进行分组，最后不足5个的忽略。对每组数进行排序（如插入排序）求取其中位数。
2. 把上一步的所有中位数移到数组的前面，对这些中位数递归调用BFPRT算法求得他们的中位数。
3. 将上一步得到的中位数作为划分的主元进行整个数组的划分。
4. 判断第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个数。