本章给出了用来分析分而治之算法复杂性的数学方法,并通过推导最小最大问题和排序问题的复杂性下限来证明分而治之算法对于求解这两种问题是最优的(因为算法的复杂性与下限一致)。
2.1 算法思想
分而治之方法与软件设计的模块化方法非常相似。为了解决一个大的问题,可以: 1) 把它分成两个或多个更小的问题; 2) 分别解决每个小问题; 3) 把各小问题的解答组合起来,即可得到原问题的解答。小问题通常与原问题相似,可以递归地使用分而治之策略来解决。
例2-1 [找出伪币] 给你一个装有1 6个硬币的袋子。1 6个硬币中有一个是伪造的,并且那个伪造的硬币比真的硬币要轻一些。你的任务是找出这个伪造的硬币。为了帮助你完成这一任务,将提供一台可用来比较两组硬币重量的仪器,利用这台仪器,可以知道两组硬币的重量是否相同。比较硬币1与硬币2的重量。假如硬币1比硬币2轻,则硬币1是伪造的;假如硬币2比硬币1轻,则硬币2是伪造的。这样就完成了任务。假如两硬币重量相等,则比较硬币3和硬币4。同样,假如有一个硬币轻一些,则寻找伪币的任务完成。假如两硬币重量相等,则继续比较硬币5和硬币6。按照这种方式,可以最多通过8次比较来判断伪币的存在并找出这一伪币。
另外一种方法就是利用分而治之方法。假如把1 6硬币的例子看成一个大的问题。第一步,把这一问题分成两个小问题。随机选择8个硬币作为第一组称为A组,剩下的8个硬币作为第二组称为B组。这样,就把1 6个硬币的问题分成两个8硬币的问题来解决。第二步,判断A和B组中是否有伪币。可以利用仪器来比较A组硬币和B组硬币的重量。假如两组硬币重量相等,则可以判断伪币不存在。假如两组硬币重量不相等,则存在伪币,并且可以判断它位于较轻的那一组硬币中。最后,在第三步中,用第二步的结果得出原先1 6个硬币问题的答案。若仅仅判断硬币是否存在,则第三步非常简单。无论A组还是B组中有伪币,都可以推断这1 6个硬币中存在伪币。因此,仅仅通过一次重量的比较,就可以判断伪币是否存在。
现在假设需要识别出这一伪币。把两个或三个硬币的情况作为不可再分的小问题。注意如果只有一个硬币,那么不能判断出它是否就是伪币。在一个小问题中,通过将一个硬币分别与其他两个硬币比较,最多比较两次就可以找到伪币。这样,1 6硬币的问题就被分为两个8硬币(A组和B组)的问题。通过比较这两组硬币的重量,可以判断伪币是否存在。如果没有伪币,则算法终止。否则,继续划分这两组硬币来寻找伪币。假设B是轻的那一组,因此再把它分成两组,每组有4个硬币。称其中一组为B1,另一组为B2。比较这两组,肯定有一组轻一些。如果B1轻,则伪币在B1中,再将B1又分成两组,每组有两个硬币,称其中一组为B1a,另一组为B1b。比较这两组,可以得到一个较轻的组。由于这个组只有两个硬币,因此不必再细分。比较组中两个硬币的重量,可以立即知道哪一个硬币轻一些。较轻的硬币就是所要找的伪币。
例2-2 [金块问题] 有一个老板有一袋金块。每个月将有两名雇员会因其优异的表现分别被奖励一个金块。按规矩,排名第一的雇员将得到袋中最重的金块,排名第二的雇员将得到袋中最轻的金块。根据这种方式,除非有新的金块加入袋中,否则第一名雇员所得到的金块总是比第二名雇员所得到的金块重。如果有新的金块周期性的加入袋中,则每个月都必须找出最轻和最重的金块。假设有一台比较重量的仪器,我们希望用最少的比较次数找出最轻和最重的金块。
假设袋中有n 个金块。可以用函数M a x(程序1 - 3 1)通过n-1次比较找到最重的金块。找到最重的金块后,可以从余下的n-1个金块中用类似的方法通过n-2次比较找出最轻的金块。这样,比较的总次数为2n-3。程序2 - 2 6和2 - 2 7是另外两种方法,前者需要进行2n-2次比较,后者最多需要进行2n-2次比较。
下面用分而治之方法对这个问题进行求解。当n很小时,比如说, n≤2,识别出最重和最轻的金块,一次比较就足够了。当n 较大时(n>2),第一步,把这袋金块平分成两个小袋A和B。第二步,分别找出在A和B中最重和最轻的金块。设A中最重和最轻的金块分别为HA 与LA,以此类推,B中最重和最轻的金块分别为HB 和LB。第三步,通过比较HA 和HB,可以找到所有金块中最重的;通过比较LA 和LB,可以找到所有金块中最轻的。在第二步中,若n>2,则递归地应用分而治之方法。
假设n= 8。这个袋子被平分为各有4个金块的两个袋子A和B。为了在A中找出最重和最轻的金块,A中的4个金块被分成两组A1和A2。每一组有两个金块,可以用一次比较在A中找出较重的金块HA1和较轻的金块LA1。经过另外一次比较,又能找出HA 2和LA 2。现在通过比较HA1和HA2,能找出HA;通过LA 1和LA2的比较找出LA。这样,通过4次比较可以找到HA 和LA。同样需要另外4次比较来确定HB 和LB。通过比较HA 和HB(LA 和LB),就能找出所有金块中最重和最轻的。因此,当n= 8时,这种分而治之的方法需要1 0次比较。如果使用程序1 - 3 1,则需要1 3次比较。如果使用程序2 - 2 6和2 - 2 7,则最多需要1 4次比较。设c(n)为使用分而治之方法所需要的比较次数。为了简便,假设n是2的幂。当n= 2时,c(n) = 1。对于较大的n,c(n) = 2c(n/ 2 ) + 2。当n是2的幂时,使用迭代方法(见例2 - 2 0)可知
c(n) = 3n/ 2 - 2。在本例中,使用分而治之方法比逐个比较的方法少用了2 5%的比较次数。
例2-3 [矩阵乘法] 两个n×n 阶的矩阵A与B的乘积是另一个n×n 阶矩阵C,C可表示为假如每一个C(i, j) 都用此公式计算,则计算C所需要的操作次数为n3 m+n2 (n- 1) a,其中m表示一次乘法,a 表示一次加法或减法。
为了得到两个矩阵相乘的分而治之算法,需要: 1) 定义一个小问题,并指明小问题是如何进行乘法运算的; 2) 确定如何把一个大的问题划分成较小的问题,并指明如何对这些较小的问题进行乘法运算; 3) 最后指出如何根据小问题的结果得到大问题的结果。为了使讨论简便,假设n 是2的幂(也就是说, n是1,2,4,8,1 6,.)。
首先,假设n= 1时是一个小问题,n> 1时为一个大问题。后面将根据需要随时修改这个假设。对于1×1阶的小矩阵,可以通过将两矩阵中的两个元素直接相乘而得到结果。
考察一个n> 1的大问题。可以将这样的矩阵分成4个n/ 2×n/ 2阶的矩阵A1,A2,A3,和A4。当n 大于1且n 是2的幂时,n/ 2也是2的幂。因此较小矩阵也满足前面对矩阵大小的假设。矩阵Bi 和Ci 的定义与此类似.
根据上述公式,经过8次n/ 2×n/ 2阶矩阵乘法和4次n/ 2×n/ 2阶矩阵的加法,就可以计算出A与B的乘积。因此,这些公式能帮助我们实现分而治之算法。在算法的第二步,将递归使用分而治之算法把8个小矩阵再细分(见程序2 - 1 9)。算法的复杂性为(n3 ),此复杂性与程序2 - 2 4直接使用公式(2 - 1)所得到的复杂性是一样的。事实上,由于矩阵分割和再组合所花费的额外开销,使用分而治之算法得出结果的时间将比用程序2 - 2 4还要长。
为了得到更快的算法,需要简化矩阵分割和再组合这两个步骤。一种方案是使用S t r a s s e n方法得到7个小矩阵。这7个小矩阵为矩阵D, E, ., J,矩阵D到J可以通过7次矩阵乘法, 6次矩阵加法,和4次矩阵减法计算得出。前述的4个小矩阵可以由矩阵D到J通过6次矩阵加法和两次矩阵减法得出.
用上述方案来解决n= 2的矩阵乘法。将某矩阵A和B相乘得结果C,如下所示:
因为n> 1,所以将A、B两矩阵分别划分为4个小矩阵,每个矩阵为1×1阶,仅包含一个元素。1×1阶矩阵的乘法为小问题,因此可以直接进行运算。利用计算D~J的公式,得:
D= 1(6-8)=-2
E= 4(7-5)= 8
F=(3 + 4)5 = 3 5
G=(1 + 2)8 = 2 4
H=(3-1)(5 + 6)= 2 2
I=(2-4)(7 + 8)=-3 0
J=(1 + 4)(5 + 8)= 6 5
根据以上结果可得:
对于上面这个2×2的例子,使用分而治之算法需要7次乘法和1 8次加/减法运算。而直接使用公式(2 - 1),则需要8次乘法和7次加/减法。要想使分而治之算法更快一些,则一次乘法所花费的时间必须比11次加/减法的时间要长。
假定S t r a s s e n矩阵分割方案仅用于n≥8的矩阵乘法,而对于n<8的矩阵乘法则直接利用公式(2 - 1)进行计算。则n= 8时,8×8矩阵相乘需要7次4×4矩阵乘法和1 8次4×4矩阵加/减法。每次矩阵乘法需花费6 4m+ 4 8a次操作,每次矩阵加法或减法需花费1 6a次操作。因此总的操作次数为7 ( 6 4m+ 4 8a) + 1 8 ( 1 6a) = 4 4 8m+ 6 2 4a。而使用直接计算方法,则需要5 1 2m+ 4 4 8a次操作。要使S t r a s s e n方法比直接计算方法快,至少要求5 1 2-4 4 8次乘法的开销比6 2 4-4 4 8次加/减法的开销大。或者说一次乘法的开销应该大于近似2 . 7 5次加/减法的开销。
假定n<1 6的矩阵是一个“小”问题,S t r a s s e n的分解方案仅仅用于n≥1 6的情况,对于n<1 6的矩阵相乘,直接利用公式( 2 - 1)。则当n= 1 6时使用分而治之算法需要7 ( 5 1 2m+ 4 4 8a) +1 8 ( 6 4a) = 3 5 8 4m+ 4 2 8 8a次操作。直接计算时需要4 0 9 6m+ 3 8 4 0a次操作。若一次乘法的开销与一次加/减法的开销相同,则S t r a s s e n方法需要7 8 7 2次操作及用于问题分解的额外时间,而直接计算方法则需要7 9 3 6次操作加上程序中执行f o r循环以及其他语句所花费的时间。即使直接计算方法所需要的操作次数比St r a s s e n方法少,但由于直接计算方法需要更多的额外开销,因此它也不见得会比S t r a s s e n方法快。
n 的值越大,Strassen 方法与直接计算方法所用的操作次数的差异就越大,因此对于足够大的n,Strassen 方法将更快。设t (n) 表示使用Strassen 分而治之方法所需的时间。因为大的矩阵会被递归地分割成小矩阵直到每个矩阵的大小小于或等于k(k至少为8,也许更大,具体值由计算机的性能决定). 用迭代方法计算,可得t(n) = (nl og27 )。因为l og27 ≈2 . 8 1,所以与直接计算方法的复杂性(n3 )相比,分而治之矩阵乘法算法有较大的改进。
注意事项
分而治之方法很自然地导致了递归算法的使用。在许多例子里,这些递归算法在递归程序中得到了很好的运用。实际上,在许多情况下,所有为了得到一个非递归程序的企图都会导致采用一个模拟递归栈。不过在有些情况下,不使用这样的递归栈而采用一个非递归程序来完成分而治之算法也是可能的,并且在这种方式下,程序得到结果的速度会比递归方式更快。解决金块问题的分而治之算法(例2 - 2)和归并排序方法( 2 . 3节)就可以不利用递归而通过一个非递归程序来更快地完成。
例2-4 [金块问题] 用例2 - 2的算法寻找8个金块中最轻和最重金块的工作可以用二叉树来表示。这棵树的叶子分别表示8个金块(a, b,., h),每个阴影节点表示一个包含其子树中所有叶子的问题。因此,根节点A表示寻找8个金块中最轻、最重金块的问题,而节点B表示找出a,b,c 和d 这4个金块中最轻和最重金块的问题。算法从根节点开始。由根节点表示的8金块问题被划分成由节点B和C所表示的两个4金块问题。在B节点,4金块问题被划分成由D和E所表示的2金块问题。可通过比较金块a 和b 哪一个较重来解决D节点所表示的2金块问题。在解决了D和E所表示的问题之后,可以通过比较D和E中所找到的轻金块和重金块来解决B表示的问题。接着在F,G和C上重复这一过程,最后解决问题A。
可以将递归的分而治之算法划分成以下的步骤:
1) 从图2 - 2中的二叉树由根至叶的过程中把一个大问题划分成许多个小问题,小问题的大小为1或2。
2) 比较每个大小为2的问题中的金块,确定哪一个较重和哪一个较轻。在节点D、E、F和G上完成这种比较。大小为1的问题中只有一个金块,它既是最轻的金块也是最重的金块。
3) 对较轻的金块进行比较以确定哪一个金块最轻,对较重的金块进行比较以确定哪一个金块最重。对于节点A到C执行这种比较。
根据上述步骤,可以得出程序1 4 - 1的非递归代码。该程序用于寻找到数组w [ 0 : n - 1 ]中的最小数和最大数,若n < 1,则程序返回f a l s e,否则返回t r u e。
当n≥1时,程序1 4 - 1给M i n和M a x置初值以使w [ M i n ]是最小的重量,w [ M a x ]为最大的重量。
首先处理n≤1的情况。若n>1且为奇数,第一个重量w [ 0 ]将成为最小值和最大值的候选值,因此将有偶数个重量值w [ 1 : n - 1 ]参与f o r循环。当n 是偶数时,首先将两个重量值放在for 循环外进行比较,较小和较大的重量值分别置为Min和Max,因此也有偶数个重量值w[2:n-1]参与for循环。
在for 循环中,外层if 通过比较确定( w [ i ] , w [ i + 1 ] )中的较大和较小者。此工作与前面提到的分而治之算法步骤中的2) 相对应,而内层的i f负责找出较小重量值和较大重量值中的最小值和
最大值,这个工作对应于3 )。for 循环将每一对重量值中较小值和较大值分别与当前的最小值w [ M i n ]和最大值w [ M a x ]进行比较,根据比较结果来修改M i n和M a x(如果必要)。
下面进行复杂性分析。注意到当n为偶数时,在for 循环外部将执行一次比较而在f o r循环内部执行3 ( n / 2 - 1 )次比较,比较的总次数为3 n / 2 - 2。当n 为奇数时,f o r循环外部没有执行比较,而内部执行了3(n-1)/2次比较。因此无论n 为奇数或偶数,当n>0时,比较的总次数为「3n/2ù-2次。
程序14-1 找出最小值和最大值的非递归程序
template
bool MinMax(T w[], int n, T& Min, T& Max)
{// 寻找w [ 0 : n - 1 ]中的最小和最大值
// 如果少于一个元素,则返回f a l s e
// 特殊情形: n <= 1
if (n < 1) return false;
if (n == 1) {Min = Max = 0;
return true;}
/ /对Min 和M a x进行初始化
int s; // 循环起点
if (n % 2) {// n 为奇数
Min = Max = 0;
s = 1;}
else {// n为偶数,比较第一对
if (w[0] > w[1]) {
Min = 1;
Max = 0;}
else {Min = 0;
Max = 1;}
s = 2;}
// 比较余下的数对
for (int i = s; i < n; i += 2) {
// 寻找w[i] 和w [ i + 1 ]中的较大者
// 然后将较大者与w [ M a x ]进行比较
// 将较小者与w [ M i n ]进行比较
if (w[i] > w[i+1]) {
if (w[i] > w[Max]) Max = i;
if (w[i+1] < w[Min]) Min = i + 1;}
else {
if (w[i+1] > w[Max]) Max = i + 1;
if (w[i] < w[Min]) Min = i;}
}
return true;
}
练习
1. 将例1 4 - 1的分而治之算法扩充到n> 1个硬币的情形。需要进行多少次重量的比较?
2. 考虑例1 4 - 1的伪币问题。假设把条件“伪币比真币轻”改为“伪币与真币的重量不同”,同样假定袋中有n 个硬币。
1) 给出相应分而治之算法的形式化描述,该算法可输出信息“不存在伪币”或找出伪币。算法应递归地将大的问题划分成两个较小的问题。需要多少次比较才能找到伪币(如果存在伪币)?
2) 重复1) ,但把大问题划分为三个较小问题。
3. 1) 编写一个C++ 程序,实现例1 4 - 2中寻找n 个元素中最大值和最小值的两种方案。使用递归来完成分而治之方案。
2) 程序2 - 2 6和2 - 2 7是另外两个寻找n 个元素中最大值和最小值的代码。试分别计算出每段程序所需要的最少和最大比较次数。
3) 在n 分别等于1 0 0,1 0 0 0或10 000的情况下,比较1)、2)中的程序和程序1 4 - 1的运行时间。对于程序2 - 2 7,使用平均时间和最坏情况下的时间。1)中的程序和程序2 - 2 6应具有相同的平均时间和最坏情况下的时间。
4) 注意到如果比较操作的开销不是很高,分而治之算法在最坏情况下不会比其他算法优越,为什么?它的平均时间优于程序2 - 2 7吗?为什么?
4. 证明直接运用公式(1 4 -2)~(1 4 - 5)得出结果的矩阵乘法的分而治之算法的复杂性为(n3 )。因此相应的分而治之程序将比程序2 - 2 4要慢。
5. 用迭代的方法来证明公式(1 4 - 6)的递归值为(n l og27)。
*6. 编写S t r a s s e n矩阵乘法程序。利用不同的k 值(见公式(1 4 - 6))进行实验,以确定k 为何值时程序性能最佳。比较程序及程序2 - 2 4的运行时间。可取n 为2的幂来进行比较。
7. 当n 不是2的幂时,可以通过增加矩阵的行和列来得到一个大小为2的幂的矩阵。假设使用最少的行数和列数将矩阵扩充为m 阶矩阵,其中m 为2的幂。
1) 求m / n。
2) 可使用哪些矩阵项组成新的行和列,以使新矩阵A' 和B' 相乘时,原来的矩阵A和B相乘的结果会出现在C' 的左上角?
3) 使用S t r a s s e n方法计算A' * B' 所需要的时间为(m2.81 )。给出以n 为变量的运行时间表达式。
2.2 应用
2.2.1 残缺棋盘
残缺棋盘(defective chessboard)是一个有2k×2k 个方格的棋盘,其中恰有一个方格残缺。图2 - 3给出k≤2时各种可能的残缺棋盘,其中残缺的方格用阴影表示。注意当k= 0时,仅存在一种可能的残缺棋盘(如图1 4 - 3 a所示)。事实上,对于任意k,恰好存在22k 种不同的残缺棋盘。
残缺棋盘的问题要求用三格板(t r i o m i n o e s)覆盖残缺棋盘(如图1 4 - 4所示)。在此覆盖中,两个三格板不能重叠,三格板不能覆盖残缺方格,但必须覆盖其他所有的方格。在这种限制条件下,所需要的三格板总数为( 22k -1 ) / 3。可以验证( 22k -1 ) / 3是一个整数。k 为0的残缺棋盘很容易被覆盖,因为它没有非残缺的方格,用于覆盖的三格板的数目为0。当k= 1时,正好存在3个非残缺的方格,并且这三个方格可用图1 4 - 4中的某一方向的三格板来覆盖。
用分而治之方法可以很好地解决残缺棋盘问题。这一方法可将覆盖2k×2k 残缺棋盘的问题转化为覆盖较小残缺棋盘的问题。2k×2k 棋盘一个很自然的划分方法就是将它划分为如图1 4 - 5 a所示的4个2k - 1×2k - 1 棋盘。注意到当完成这种划分后, 4个小棋盘中仅仅有一个棋盘存在残缺方格(因为原来的2k×2k 棋盘仅仅有一个残缺方格)。首先覆盖其中包含残缺方格的2k - 1×2k - 1 残缺棋盘,然后把剩下的3个小棋盘转变为残缺棋盘,为此将一个三格板放在由这3个小棋盘形成的角上,如图14-5b 所示,其中原2k×2k 棋盘中的残缺方格落入左上角的2k - 1×2k - 1 棋盘。可以采用这种分割技术递归地覆盖2k×2k 残缺棋盘。当棋盘的大小减为1×1时,递归过程终止。此时1×1的棋盘中仅仅包含一个方格且此方格残缺,所以无需放置三格板。
可以将上述分而治之算法编写成一个递归的C++ 函数Ti l e B o a r d (见程序1 4 - 2 )。该函数定义了一个全局的二维整数数组变量B o a r d来表示棋盘。B o a r d [ 0 ] [ 0 ]表示棋盘中左上角的方格。该函数还定义了一个全局整数变量t i l e,其初始值为0。函数的输入参数如下:
? tr 棋盘中左上角方格所在行。
? tc 棋盘中左上角方格所在列。
? dr 残缺方块所在行。
? dl 残缺方块所在列。
? size 棋盘的行数或列数。
Ti l e B o a r d函数的调用格式为Ti l e B o a r d(0,0, dr, dc,size),其中s i z e = 2k。覆盖残缺棋盘所需要的三格板数目为( s i z e2 -1 ) / 3。函数TileBoard 用整数1到( s i z e2-1 ) / 3来表示这些三格板,并用三格板的标号来标记被该三格板覆盖的非残缺方格。
令t (k) 为函数Ti l e B o a r d覆盖一个2k×2k 残缺棋盘所需要的时间。当k= 0时,s i z e等于1,覆盖它将花费常数时间d。当k > 0时,将进行4次递归的函数调用,这些调用需花费的时间为4t (k-1 )。除了这些时间外, if 条件测试和覆盖3个非残缺方格也需要时间,假设用常数c 表示这些额外时间。可以得到以下递归表达式:
程序14-2 覆盖残缺棋盘
void TileBoard(int tr, int tc, int dr, int dc, int size)
{// 覆盖残缺棋盘
if (size == 1) return;
int t = tile++, // 所使用的三格板的数目
s = size/2; // 象限大小
/ /覆盖左上象限
if (dr < tr + s && dc < tc + s)
// 残缺方格位于本象限
Ti l e B o a r d ( t r, tc, dr, dc, s);
else {// 本象限中没有残缺方格
// 把三格板t 放在右下角
Board[tr + s - 1][tc + s - 1] = t;
// 覆盖其余部分
Ti l e B o a r d ( t r, tc, tr+s-1, tc+s-1, s);}
/ /覆盖右上象限
if (dr < tr + s && dc >= tc + s)
// 残缺方格位于本象限
Ti l e B o a r d ( t r, tc+s, dr, dc, s);
else {// 本象限中没有残缺方格
// 把三格板t 放在左下角
Board[tr + s - 1][tc + s] = t;
// 覆盖其余部分
Ti l e B o a r d ( t r, tc+s, tr+s-1, tc+s, s);}
/ /覆盖左下象限
if (dr >= tr + s && dc < tc + s)
// 残缺方格位于本象限
TileBoard(tr+s, tc, dr, dc, s);
else {// 把三格板t 放在右上角
Board[tr + s][tc + s - 1] = t;
// 覆盖其余部分
TileBoard(tr+s, tc, tr+s, tc+s-1, s);}
// 覆盖右下象限
if (dr >= tr + s && dc >= tc + s)
// 残缺方格位于本象限
TileBoard(tr+s, tc+s, dr, dc, s);
else {// 把三格板t 放在左上角
Board[tr + s][tc + s] = t;
// 覆盖其余部分
TileBoard(tr+s, tc+s, tr+s, tc+s, s);}
}
void OutputBoard(int size)
{
for (int i = 0; i < size; i++) {
for (int j = 0; j < size; j++)
cout << setw (5) << Board[i][j];
cout << endl;
}
}
可以用迭代的方法来计算这个表达式(见例2 - 2 0),可得t (k )= ( 4k )= (所需的三格板的数目)。由于必须花费至少( 1 )的时间来放置每一块三格表,因此不可能得到一个比分而治之算法更快的算法。
2.2.2 归并排序
可以运用分而治之方法来解决排序问题,该问题是将n 个元素排成非递减顺序。分而治之方法通常用以下的步骤来进行排序算法:若n 为1,算法终止;否则,将这一元素集合分割成两个或更多个子集合,对每一个子集合分别排序,然后将排好序的子集合归并为一个集合。
假设仅将n 个元素的集合分成两个子集合。现在需要确定如何进行子集合的划分。一种可能性就是把前面n- 1个元素放到第一个子集中(称为A),最后一个元素放到第二个子集里(称为B)。按照这种方式对A递归地进行排序。由于B仅含一个元素,所以它已经排序完毕,在A排完序后,只需要用程序2 - 1 0中的函数i n s e r t将A和B合并起来。把这种排序算法与I n s e r t i o n S o r t(见程序2 - 1 5)进行比较,可以发现这种排序算法实际上就是插入排序的递归算法。该算法的复杂性为O (n 2 )。把n 个元素划分成两个子集合的另一种方法是将含有最大值的元素放入B,剩下的放入A中。然后A被递归排序。为了合并排序后的A和B,只需要将B添加到A中即可。假如用函数M a x(见程序1 - 3 1)来找出最大元素,这种排序算法实际上就是S e l e c t i o n S o r t(见程序2 - 7)的递归算法。
假如用冒泡过程(见程序2 - 8)来寻找最大元素并把它移到最右边的位置,这种排序算法就是B u b b l e S o r t(见程序2 - 9)的递归算法。这两种递归排序算法的复杂性均为(n2 )。若一旦发现A已经被排好序就终止对A进行递归分割,则算法的复杂性为O(n2 )(见例2 - 1 6和2 - 1 7)。
上述分割方案将n 个元素分成两个极不平衡的集合A和B。A有n- 1个元素,而B仅含一个元素。下面来看一看采用平衡分割法会发生什么情况: A集合中含有n/k 个元素,B中包含其余的元素。递归地使用分而治之方法对A和B进行排序。然后采用一个被称之为归并( m e rg e)的过程,将已排好序的A和B合并成一个集合。
例2-5 考虑8个元素,值分别为[ 1 0,4,6,3,8,2,5,7 ]。如果选定k = 2,则[ 1 0 , 4 , 6 , 3 ]和[ 8 , 2 , 5 , 7 ]将被分别独立地排序。结果分别为[ 3 , 4 , 6 , 1 0 ]和[ 2 , 5 , 7 , 8 ]。从两个序列的头部开始归并这两个已排序的序列。元素2比3更小,被移到结果序列;3与5进行比较,3被移入结果序列;4与5比较,4被放入结果序列;5和6比较,.。如果选择k= 4,则序列[ 1 0 , 4 ]和[ 6 , 3 , 8 , 2 , 5 , 7 ]将被排序。排序结果分别为[ 4 , 1 0 ]和[ 2 , 3 , 5 , 6 , 7 , 8 ]。当这两个排好序的序列被归并后,即可得所需要的排序序列。
图2 - 6给出了分而治之排序算法的伪代码。算法中子集合的数目为2,A中含有n/k个元素。
template
void sort( T E, int n)
{ / /对E中的n 个元素进行排序, k为全局变量
if (n >= k) {
i = n/k;
j = n-i;
令A 包含E中的前i 个元素
令B 包含E中余下的j 个元素
s o r t ( A , i ) ;
s o r t ( B , j ) ;
m e rge(A,B,E,i,j,); //把A 和B 合并到E
}
else 使用插入排序算法对E 进行排序
}
图14-6 分而治之排序算法的伪代码
从对归并过程的简略描述中,可以明显地看出归并n个元素所需要的时间为O (n)。设t (n)为分而治之排序算法(如图1 4 - 6所示)在最坏情况下所需花费的时间,则有以下递推公式:
其中c 和d 为常数。当n / k≈n-n / k 时,t (n) 的值最小。因此当k= 2时,也就是说,当两个子集合所包含的元素个数近似相等时, t (n) 最小,即当所划分的子集合大小接近时,分而治之算法通常具有最佳性能。
可以用迭代方法来计算这一递推方式,结果为t(n)= (nl o gn)。虽然这个结果是在n为2的幂时得到的,但对于所有的n,这一结果也是有效的,因为t(n) 是n 的非递减函数。t(n) =(nl o gn) 给出了归并排序的最好和最坏情况下的复杂性。由于最好和最坏情况下的复杂性是一样的,因此归并排序的平均复杂性为t (n)= (nl o gn)。
图2 - 6中k= 2的排序方法被称为归并排序( m e rge sort ),或更精确地说是二路归并排序(two-way merge sort)。下面根据图1 4 - 6中k= 2的情况(归并排序)来编写对n 个元素进行排序的C + +函数。一种最简单的方法就是将元素存储在链表中(即作为类c h a i n的成员(程序3 -8))。在这种情况下,通过移到第n/ 2个节点并打断此链,可将E分成两个大致相等的链表。
归并过程应能将两个已排序的链表归并在一起。如果希望把所得到C + +程序与堆排序和插入排序进行性能比较,那么就不能使用链表来实现归并排序,因为后两种排序方法中都没有使用链表。为了能与前面讨论过的排序函数作比较,归并排序函数必须用一个数组a来存储元素集合E,并在a 中返回排序后的元素序列。为此按照下述过程来对图1 4 - 6的伪代码进行细化:当集合E被化分成两个子集合时,可以不必把两个子集合的元素分别复制到A和B中,只需简单地在集合E中保持两个子集合的左右边界即可。接下来对a 中的初始序列进行排序,并将所得到的排序序列归并到一个新数组b中,最后将它们复制到a 中。图1 4 - 6的改进版见图1 4 - 7。
template
M e rgeSort( T a[], int left, int right)
{ / /对a [ l e f t : r i g h t ]中的元素进行排序
if (left < right) {//至少两个元素
int i = (left + right)/2; //中心位置
M e rgeSort(a, left, i);
M e rgeSort(a, i+1, right);
M e rge(a, b, left, i, right); //从a 合并到b
Copy(b, a, left, right); //结果放回a
}
}
图14-7 分而治之排序算法的改进
可以从很多方面来改进图1 4 - 7的性能,例如,可以容易地消除递归。如果仔细地检查图1 4 - 7中的程序,就会发现其中的递归只是简单地重复分割元素序列,直到序列的长度变成1为止。当序列的长度变为1时即可进行归并操作,这个过程可以用n 为2的幂来很好地描述。长度为1的序列被归并为长度为2的有序序列;长度为2的序列接着被归并为长度为4的有序序列;这个过程不断地重复直到归并为长度为n 的序列。图1 4 - 8给出n= 8时的归并(和复制)过程,方括号表示一个已排序序列的首和尾。
初始序列[8] [4] [5] [6] [2] [1] [7] [3]
归并到b [4 8] [5 6] [1 2] [3 7]
复制到a [4 8] [5 6] [1 2] [3 7]
归并到b [4 5 6 8] [1 2 3 7]
复制到a [4 5 6 8] [1 2 3 7]
归并到b [1 2 3 4 5 6 7 8]
复制到a [1 2 3 4 5 6 7 8]
图14-8 归并排序的例子
另一种二路归并排序算法是这样的:首先将每两个相邻的大小为1的子序列归并,然后对上一次归并所得到的大小为2的子序列进行相邻归并,如此反复,直至最后归并到一个序列,归并过程完成。通过轮流地将元素从a 归并到b 并从b 归并到a,可以虚拟地消除复制过程。二路归并排序算法见程序1 4 - 3。
程序14-3 二路归并排序
template
void MergeSort(T a[], int n)
{// 使用归并排序算法对a[0:n-1] 进行排序
T *b = new T [n];
int s = 1; // 段的大小
while (s < n) {
MergePass(a, b, s, n); // 从a归并到b
s += s;
MergePass(b, a, s, n); // 从b 归并到a
s += s;
}
}
为了完成排序代码,首先需要完成函数M e rg e P a s s。函数M e rg e P a s s(见程序1 4 - 4)仅用来确定欲归并子序列的左端和右端,实际的归并工作由函数M e rg e (见程序1 4 - 5 )来完成。函数M e rg e要求针对类型T定义一个操作符< =。如果需要排序的数据类型是用户自定义类型,则必须重载操作符< =。这种设计方法允许我们按元素的任一个域进行排序。重载操作符< =的目的是用来比较需要排序的域。
程序14-4 MergePass函数
template
void MergePass(T x[], T y[], int s, int n)
{// 归并大小为s的相邻段
int i = 0;
while (i <= n - 2 * s) {
// 归并两个大小为s的相邻段
Merge(x, y, i, i+s-1, i+2*s-1);
i = i + 2 * s;
}
// 剩下不足2个元素
if (i + s < n) Merge(x, y, i, i+s-1, n-1);
else for (int j = i; j <= n-1; j++)
// 把最后一段复制到y
y[j] = x[j];
}
程序14-5 Merge函数
template
void Merge(T c[], T d[], int l, int m, int r)
{// 把c[l:m]] 和c[m:r] 归并到d [ l : r ] .
int i = l, // 第一段的游标
j = m+1, // 第二段的游标
k = l; // 结果的游标
/ /只要在段中存在i和j,则不断进行归并
while ((i <= m) && (j <= r))
if (c[i] <= c[j]) d[k++] = c[i++];
else d[k++] = c[j++];
// 考虑余下的部分
if (i > m) for (int q = j; q <= r; q++)
d[k++] = c[q];
else for (int q = i; q <= m; q++)
d[k++] = c[q];
}
自然归并排序(natural merge sort)是基本归并排序(见程序1 4 - 3)的一种变化。它首先对输入序列中已经存在的有序子序列进行归并。例如,元素序列[ 4,8,3,7,1,5,6,2 ]中包含有序的子序列[ 4,8 ],[ 3,7 ],[ 1,5,6 ]和[ 2 ],这些子序列是按从左至右的顺序对元素表进行扫描而产生的,若位置i 的元素比位置i+ 1的元素大,则从位置i 进行分割。对于上面这个元素序列,可找到四个子序列,子序列1和子序列2归并可得[ 3 , 4 , 7 , 8 ],子序列3和子序列4归并可得[ 1 , 2 , 5 , 6 ],最后,归并这两个子序列得到[ 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 ]。因此,对于上述元素序列,仅仅使用了两趟归并,而程序1 4 - 3从大小为1的子序列开始,需使用三趟归并。作为一个极端的例子,假设输入的元素序列已经排好序并有n个元素,自然归并排序法将准确地识别该序列不必进行归并排序,但程序1 4 - 3仍需要进行[ l o g2 n] 趟归并。因此自然归并排序将在(n) 的时间内完成排序。而程序1 4 - 3将花费(n l o gn) 的时间。
2.2.3 快速排序
分而治之方法还可以用于实现另一种完全不同的排序方法,这种排序法称为快速排序(quick sort)。在这种方法中, n 个元素被分成三段(组):左段l e f t,右段r i g h t和中段m i d d l e。中段仅包含一个元素。左段中各元素都小于等于中段元素,右段中各元素都大于等于中段元素。因此l e f t和r i g h t中的元素可以独立排序,并且不必对l e f t和r i g h t的排序结果进行合并。m i d d l e中的元素被称为支点( p i v o t )。图1 4 - 9中给出了快速排序的伪代码。
/ /使用快速排序方法对a[ 0 :n- 1 ]排序
从a[ 0 :n- 1 ]中选择一个元素作为m i d d l e,该元素为支点
把余下的元素分割为两段left 和r i g h t,使得l e f t中的元素都小于等于支点,而right 中的元素都大于等于支点
递归地使用快速排序方法对left 进行排序
递归地使用快速排序方法对right 进行排序
所得结果为l e f t + m i d d l e + r i g h t
图14-9 快速排序的伪代码
考察元素序列[ 4 , 8 , 3 , 7 , 1 , 5 , 6 , 2 ]。假设选择元素6作为支点,则6位于m i d d l e;4,3,1,5,2位于l e f t;8,7位于r i g h t。当left 排好序后,所得结果为1,2,3,4,5;当r i g h t排好序后,所得结果为7,8。把right 中的元素放在支点元素之后, l e f t中的元素放在支点元素之前,即可得到最终的结果[ 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 ]。
把元素序列划分为l e f t、m i d d l e和r i g h t可以就地进行(见程序1 4 - 6)。在程序1 4 - 6中,支点总是取位置1中的元素。也可以采用其他选择方式来提高排序性能,本章稍后部分将给出这样一种选择。
程序14-6 快速排序
template
void QuickSort(T*a, int n)
{// 对a[0:n-1] 进行快速排序
{// 要求a[n] 必需有最大关键值
quickSort(a, 0, n-1);
template
void quickSort(T a[], int l, int r)
{// 排序a [ l : r ], a[r+1] 有大值
if (l >= r) return;
int i = l, // 从左至右的游标
j = r + 1; // 从右到左的游标
T pivot = a[l];
// 把左侧>= pivot的元素与右侧<= pivot 的元素进行交换
while (true) {
do {// 在左侧寻找>= pivot 的元素
i = i + 1;
} while (a[i] < pivot);
do {// 在右侧寻找<= pivot 的元素
j = j - 1;
} while (a[j] > pivot);
if (i >= j) break; // 未发现交换对象
Swap(a[i], a[j]);
}
// 设置p i v o t
a[l] = a[j];
a[j] = pivot;
quickSort(a, l, j-1); // 对左段排序
quickSort(a, j+1, r); // 对右段排序
}
若把程序1 4 - 6中d o - w h i l e条件内的<号和>号分别修改为< =和> =,程序1 4 - 6仍然正确。实验结果表明使用程序1 4 - 6的快速排序代码可以得到比较好的平均性能。为了消除程序中的递归,必须引入堆栈。不过,消除最后一个递归调用不须使用堆栈。消除递归调用的工作留作练习(练习1 3)。程序1 4 - 6所需要的递归栈空间为O (n)。若使用堆栈来模拟递归,则可以把这个空间减少为O ( l o gn)。在模拟过程中,首先对left 和right 中较小者进行排序,把较大者的边界放入堆栈中。在最坏情况下l e f t总是为空,快速排序所需的计算时间为(n2 )。在最好情况下, l e f t和r i g h t中的元素数目大致相同,快速排序的复杂性为(nl o gn)。令人吃惊的是,快速排序的平均复杂性也是(nl o gn)。
定理2-1 快速排序的平均复杂性为(nl o gn)。
证明用t (n) 代表对含有n 个元素的数组进行排序的平均时间。当n≤1时,t (n)≤d,d为某一常数。当n <1时,用s 表示左段所含元素的个数。由于在中段中有一个支点元素,因此右段中元素的个数为n-s- 1。所以左段和右段的平均排序时间分别为t (s), t (n-s- 1 )。分割数组中元素所需要的时间用cn 表示,其中c 是一个常数。因为s 有同等机会取0 ~n- 1中的任何一个值.
如对(2 - 8)式中的n 使用归纳法,可得到t (n)≤kn l o ge n,其中n> 1且k=2(c+d),e~2 . 7 1 8为自然对数的基底。在归纳开始时首先验证n= 2时公式的正确性。根据公式( 1 4 - 8),可以得到t( 2 )≤2c+ 2d≤k nl o ge 2。在归纳假设部分,假定t(n)≤kn l o ge n(当2≤n<m 时,m 是任意一个比2大的整数=.
图1 4 - 1 0对本书中所讨论的算法在平均条件下和最坏条件下的复杂性进行了比较。
方法最坏复杂性平均复杂性
冒泡排序n2 n2
计数排序n2 n2
插入排序n2 n2
选择排序n2 n2
堆排序nl o gn nl o gn
归并排序nl o gn nl o gn
快速排序n2 nl o gn
图14-10 各种排序算法的比较
中值快速排序( median-of-three quick sort)是程序1 4 - 6的一种变化,这种算法有更好的平均性能。注意到在程序1 4 - 6中总是选择a [ 1 ]做为支点,而在这种快速排序算法中,可以不必使用a [ 1 ]做为支点,而是取 中大小居中的那个元素作为支点。例如,假如有三个元素,大小分别为5,9,7,那么取7为支点。为了实现中值快速排序算法,一种最简单的方式就是首先选出中值元素并与a[1] 进行交换,然后利用程序1 4 - 6完成排序。如果a [ r ]是被选出的中值元素,那么将a[1] 与a[r] 进行交换,然后将a [ 1 ](即原来的a [ r ])赋值给程序1 4 - 6中的变量p i v o t,之后继续执行程序1 4 - 6中的其余代码。
图2 - 11中分别给出了根据实验所得到的归并排序、堆排序、插入排序、快速排序的平均时间。对于每一个不同的n, 都随机产生了至少1 0 0组整数。随机整数的产生是通过反复调用s t d l i b . h库中的r a n d o m函数来实现的。如果对一组整数进行排序的时间少于1 0个时钟滴答,则继续对其他组整数进行排序,直到所用的时间不低于1 0个时钟滴答。在图2 - 11中的数据包含产生随机整数的时间。对于每一个n,在各种排序法中用于产生随机整数及其他开销的时间是相同的。因此,图2 - 11中的数据对于比较各种排序算法是很有用的。
对于足够大的n,快速排序算法要比其他算法效率更高。从图中可以看到快速排序曲线与插入排序曲线的交点横坐标比2 0略小,可通过实验来确定这个交点横坐标的精确值。可以分别用n = 1 5 , 1 6 , 1 7 , 1 8 , 1 9进行实验,以寻找精确的交点。令精确的交点横坐标为nBr e a k。当n≤nBreak 时,插入排序的平均性能最佳。当n>nBreak 时,快速排序性能最佳。当n>nBreak 时,把插入排序与快速排序组合为一个排序函数,可以提高快速排序的性能,实现方法是把程序1 4 - 6中的以下语句:
if(l >= r)r e t u r n ;
替换为
if (r-1
这里I n s e r t i o n S o r t ( a , l , r )用来对a [ 1 : r ]进行插入排序。测量修改后的快速排序算法的性能留作练习(练习2 0)。用更小的值替换n B r e a k有可能使性能进一步提高(见练习2 0)。
大多数实验表明,当n>c时(c为某一常数),在最坏情况下归并排序的性能也是最佳的。而当n≤c时,在最坏情况下插入排序的性能最佳。通过将插入排序与归并排序混合使用,可以提高归并排序的性能(练习2 1)。
2.2.4 选择
对于给定的n 个元素的数组a [ 0 : n - 1 ],要求从中找出第k小的元素。当a [ 0 : n - 1 ]被排序时,该元素就是a [ k - 1 ]。假设n = 8,每个元素有两个域k e y和I D,其中k e y是一个整数,I D是一个字符。假设这8个元素为[ ( 1 2 ,a),( 4 ,b),( 5 ,c),( 4 ,d),( 5 ,e),( 1 0 ,f),( 2 ,g),( 2 0 ,h)], 排序后得到数组[ ( 2 ,g),( 4 ,d),( 4 ,b),( 5 ,c),( 5 ,e),( 1 0 ,f),( 1 2 ,a),( 2 0 ,h) ]。如果k = 1,返回I D为g 的元素;如果k = 8,返回I D为h 的元素;如果k = 6,返回是I D为f 的元素;如果k = 2,返回I D为d 的元素。实际上,对最后一种情况,所得到的结果可能不唯一,因为排序过程中既可能将I D为d 的元素排在a [ 1 ],也可能将I D为b 的元素排在a [ 1 ],原因是它们具有相同大小的k e y,因而两个元素中的任何一个都有可能被返回。但是无论如何,如果一个元素在k = 2时被返回,另一个就必须在k = 3时被返回。
选择问题的一个应用就是寻找中值元素,此时k = [n / 2 ]。中值是一个很有用的统计量,例如中间工资,中间年龄,中间重量。其他k值也是有用的。例如,通过寻找第n / 4 , n / 2和3 n / 4这三个元素,可将人口划分为4份。
选择问题可在O ( n l o g n )时间内解决,方法是首先对这n个元素进行排序(如使用堆排序式或归并排序),然后取出a [ k - 1 ]中的元素。若使用快速排序(如图1 4 - 11所示),可以获得更好的平均性能,尽管该算法有一个比较差的渐近复杂性O( n2 )。
可以通过修写程序1 4 - 6来解决选择问题。如果在执行两个w h i l e循环后支点元素a [ l ]被交换到a [ j ] ,那么a [ l ]是a [ l : j ]中的第j - l + 1个元素。如果要寻找的第k 个元素在a [ l : r ]中,并且j - l + 1等于k,则答案就是a [ l ];如果j - l + 1 < k,那么寻找的元素是r i g h t中的第k - j + l - 1个元素,否则要寻找的元素是left 中的第k个元素。因此,只需进行0次或1次递归调用。新代码见程序1 4 - 7。S e l e c t中的递归调用可用f o r或w h i l e循环来替代(练习2 5)。
程序14-7 寻找第k 个元素
template
T Select(T a[], int n, int k)
{// 返回a [ 0 : n - 1 ]中第k小的元素
// 假定a[n] 是一个伪最大元素
if (k < 1 || k > n) throw OutOfBounds();
return select(a, 0, n-1, k);
}
template
T select(T a[], int l, int r, int k)
{// 在a [ l : r ]中选择第k小的元素
if (l >= r) return a[l];
int i = l, // 从左至右的游标
j = r + 1; // 从右到左的游标
T pivot = a[l];
// 把左侧>= pivot的元素与右侧<= pivot 的元素进行交换
while (true) {
do {// 在左侧寻找>= pivot 的元素
i = i + 1;
} while (a[i] < pivot);
do {// 在右侧寻找<= pivot 的元素
j = j - 1;
} while (a[j] > pivot);
if (i >= j) break; // 未发现交换对象
Swap(a[i], a[j]);
}
if (j - l + 1 == k) return pivot;
// 设置p i v o t
a[l] = a[j];
a[j] = pivot;
// 对一个段进行递归调用
if (j - l + 1 < k)
return select(a, j+1, r, k-j+l-1);
else return select(a, l, j-1, k);
}
程序1 4 - 7在最坏情况下的复杂性是( n2 ),此时left 总是为空,而且第k个元素总是位于r i g h t.
如果假定n 是2的幂,则可以取消公式(2 - 1 0)中的向下取整操作符。通过使用迭代方法,可以得到t (n) = (n)。若仔细地选择支点元素,则最坏情况下的时间开销也可以变成(n)。一种选择支点元素的方法是使用“中间的中间( m e d i a n - o f - m e d i a n)”规则,该规则首先将数组a中的n 个元素分成n/r 组,r 为某一整常数,除了最后一组外,每组都有r 个元素。然后通过在每组中对r 个元素进行排序来寻找每组中位于中间位置的元素。最后根据所得到的n/r 个中间元素,递归使用选择算法,求得所需要的支点元素。
例2-6 [中间的中间] 考察如下情形:r=5, n=27, 并且a= [ 2,6,8,1,4,1 0,2 0,6,2 2,11,9,8,4,3,7,8,1 6,11,1 0,8,2,1 4,1 5,1,1 2,5,4 ]。这2 7个元素可以被分为6组[ 2 , 6 , 8 , 1 , 4 ],[ 1 0 , 2 0 , 6 , 2 2 , 11 ],[ 9 , 8 , 4 , 3 , 7 ],[ 8 , 1 6 , 11 , 1 0 , 8 ],[ 2 , 1 4 , 1 5 , 1 , 1 2 ]和[ 5 , 4 ],每组的中间元素分别为4 , 11 , 7 , 1 0 , 1 2和4。[ 4 , 11 , 7 , 1 0 , 1 2 , 4 ]的中间元素为7。这个中间元素7被取为支点元素。由此可以得到l e ft= [ 2 , 6 , 1 , 4 , 6 , 4 , 3 , 2 , 1 , 5 , 4 ],m i d d l e= [ 7 ] ,r i g h t= [ 8 , 1 0 , 2 0 , 2 2 , 11 , 9 , 8 , 8 , 1 6 , 11 , 1 0 , 8 , 1 4 , 1 5 , 1 2 ]。
如果要寻找第k个元素且k< 1 2,则仅仅需要在l e f t中寻找;如果k= 1 2,则要找的元素就是支点元素;如果k> 1 2,则需要检查r i g h t中的1 5个元素。在最后一种情况下,需在r i g h t中寻找第(k- 1 2 )个元素。
定理2-2 当按“中间的中间”规则选取支点元素时,以下结论为真:
1) 若r=9, 那么当n≥9 0时,有m a x { |l e f e|, |r i g h t| }≤7n / 8。
2) 若r= 5,且a 中所有元素都不同,那么当n≥2 4时,有max{| left |, | right | }≤3n/ 4。
证明这个定理的证明留作练习2 3。
根据定理2 - 2和程序1 4 - 7可知,如果采用“中间的中间”规则并取r= 9,则用于寻找第k个元素的时间t (n)可按如下递归公式来计算:
在上述递归公式中,假设当n<9 0时使用复杂性为nl o gn的求解算法,当n≥9 0时,采用“中间的中间”规则进行分而治之求解。利用归纳法可以证明,当n≥1时有t (n)≤7 2cn (练习2 4 )。
当元素互不相同时,可以使用r= 5来得到线性时间性能。
2.2.5 距离最近的点对
给定n 个点(xi,yi)(1≤i≤n),要求找出其中距离最近的两个点。
例14-7 假设在一片金属上钻n 个大小一样的洞,如果洞太近,金属可能会断。若知道任意两个洞的最小距离,可估计金属断裂的概率。这种最小距离问题实际上也就是距离最近的点对问题。
通过检查所有的n(n- 1 ) / 2对点,并计算每一对点的距离,可以找出距离最近的一对点。这种方法所需要的时间为(n2 )。我们称这种方法为直接方法。图1 4 - 1 3中给出了分而治之求解算法的伪代码。该算法对于小的问题采用直接方法求解,而对于大的问题则首先把它划分为两个较小的问题,其中一个问题(称为A)的大小为「n /2ù,另一个问题(称为B)的大小为「n /2ù。初始时,最近的点对可能属于如下三种情形之一: 1) 两点都在A中(即最近的点对落在A中);2) 两点都在B中;3) 一点在A,一点在B。假定根据这三种情况来确定最近点对,则最近点对是所有三种情况中距离最小的一对点。在第一种情况下可对A进行递归求解,而在第二种情况下可对B进行递归求解。
if (n较小) {用直接法寻找最近点对
R e t u r n ; }
// n较大
将点集分成大致相等的两个部分A和B
确定A和B中的最近点对
确定一点在A中、另一点在B中的最近点对
从上面得到的三对点中,找出距离最小的一对点
图14-13 寻找最近的点对
为了确定第三种情况下的最近点对,需要采用一种不同的方法。这种方法取决于点集是如何被划分成A、B的。一个合理的划分方法是从xi(中间值)处划一条垂线,线左边的点属于A,线右边的点属于B。位于垂线上的点可在A和B之间分配,以便满足A、B的大小。
例2-8 考察图14-14a 中从a到n的1 4个点。这些点标绘在图14-14b 中。中点xi = 1,垂线x = 1如图14-14b 中的虚线所示。虚线左边的点(如b, c, h, n, i)属于A,右边的点(如a, e, f, j, k, l) 属于B。d, g, m 落在垂线上,可将其中两个加入A, 另一个加入B,以便A、B中包含相同的点数。假设d ,m加入A,g加入B。
设是i 的最近点对和B的最近点对中距离较小的一对点。若第三种情况下的最近点对比小。则每一个点距垂线的距离必小于,这样,就可以淘汰那些距垂线距离≥ 的点。图1 4 - 1 5中的虚线是分割线。阴影部分以分割线为中线,宽为2 。边界线及其以外的点均被淘汰掉,只有阴影中的点被保留下来,以便确定是否存在第三类点对(对应于第三种情况)其距离小于。用RA、RB 分别表示A和B中剩下的点。如果存在点对(p,q),p?A, q?B且p, q 的距离小于,则p?RA,q?RB。可以通过每次检查RA 中一个点来寻找这样的点对。假设考察RA 中的p 点,p的y 坐标为p.y,那么只需检查RB 中满足p.y- <q.y<p.y+ 的q 点,看是否存在与p 间距小于的点。在图14-16a 中给出了包含这种q 点的RB 的范围。因此,只需将RB 中位于×2 阴影内的点逐个与p 配对,以判断p 是否是距离小于的第三类点。这个×2 区域被称为是p 的比较区(comparing region)。
例2-9 考察例2 - 8中的1 4个点。A中的最近点对为(b,h),其距离约为0 . 3 1 6。B中最近点对为(f, j),其距离为0 . 3,因此= 0 . 3。当考察是否存在第三类点时,除d, g, i, l, m 以外的点均被淘汰,因为它们距分割线x= 1的距离≥ 。RA ={d, i, m},RB= {g, l},由于d 和m 的比较区中没有点,只需考察i即可。i 的比较区中仅含点l。计算i 和l的距离,发现它小于,因此(i, l) 是最近的点对。
为了确定一个距离更小的第三类点,RA 中的每个点最多只需和RB 中的6个点比较,如图1 4 - 1 6所示。
1. 选择数据结构
为了实现图1 4 - 1 3的分而治之算法,需要确定什么是“小问题”以及如何表示点。由于集合中少于两点时不存在最近点对,因此必须保证分解过程不会产生少于两点的点集。如果将少于四点的点集做为“小问题”,就可以避免产生少于两点的点集。
每个点可有三个参数:标号, x 坐标,y 坐标。假设标号为整数,每个点可用P o i n t l类(见程序1 4 - 8)来表示。为了便于按x 坐标对各个点排序,可重载操作符<=。归并排序程序如1 4 -3所示。
程序14-8 点类
class Point1 {
friend float dist(const Point1&, const Point1&);
friend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&, Point1&, float&);
friend bool closest(Point1 *, int, Point1&, Point1&,float&);
friend void main();
p u b l i c :
int operator<=(Point1 a) const
{return (x <= a.x);}
p r i v a t e :
int ID; // 点的编号
float x, y; // 点坐标
} ;
class Point2 {
friend float dist(const Point2&, const Point2&);
friend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&, Point1&, float&);
friend bool closest(Point1 *, int, Point1&, Point1&, float&);
friend void main();
p u b l i c :
int operator<=(Point2 a) const
{return (y <= a.y);}
p r i v a t e :
int p; // 数组X中相同点的索引
float x, y; // 点坐标
} ;
所输入的n 个点可以用数组X来表示。假设X中的点已按照x 坐标排序,在分割过程中如果当前考察的点是X [l :r],那么首先计算m= (l+r) / 2,X[ l:m]中的点属于A,剩下的点属于B。计算出A和B中的最近点对之后,还需要计算RA 和RB,然后确定是否存在更近的点对,其中一点属于RA,另一点属于RB。如果点已按y 坐标排序,那么可以用一种很简单的方式来测试图1 4 - 1 6。按y 坐标排序的点保存在另一个使用类P o i n t 2 (见程序14-8) 的数组中。注意到在P o i n t 2类中,为了便于y 坐标排序,已重载了操作符<=。成员p 用于指向X中的对应点。
确定了必要的数据结构之后,再来看看所要产生的代码。首先定义一个模板函数d i s t (见程序1 4 - 9 )来计算点a, b 之间的距离。T可能是P o i n t 1或P o i n t 2,因此d i s t必须是P o i n t 1和P o i n t 2类的友元。
程序14-9 计算两点距离
template
inline float dist(const T& u, const T& v)
{ / /计算点u 和v之间的距离
float dx = u.x-v. x ;
float dy = u.y-v. y ;
return sqrt(dx * dx + dy * dy);
}
如果点的数目少于两个,则函数c l o s e s t (见程序1 4 - 1 0 )返回f a l s e,如果成功时函数返回t r u e。当函数成功时,在参数a 和b 中返回距离最近的两个点,在参数d 中返回距离。代码首先验证至少存在两点,然后使用M e rg e S o r t函数(见程序14-3) 按x 坐标对X中的点排序。接下来把这些点复制到数组Y中并按y 坐标进行排序。排序完成时,对任一个i,有Y [i ] . y≤Y [i+ 1 ] . y,并且Y [i ] .p给出了点i 在X中的位置。上述准备工作做完以后,调用函数close (见程序1 4 - 11 ),该函数实际求解最近点对。
程序14-10 预处理及调用c l o s e
bool closest(Point1 X[], int n, Point1& a, Point1& b, float& d)
{// 在n >= 2 个点中寻找最近点对
// 如果少于2个点,则返回f a l s e
// 否则,在a 和b中返回距离最近的两个点
if (n < 2) return false;
// 按x坐标排序
M e r g e S o r t ( X , n ) ;
// 创建一个按y坐标排序的点数组
Point2 *Y = new Point2 [n];
for (int i = 0; i < n; i++) {
// 将点i 从X 复制到Y
Y[i].p = i;
Y[i].x = X[i].x;
Y[i].y = X[i].y;
}
M e r g e S o r t ( Y,n); // 按y坐标排序
// 创建临时数组
Point2 *Z = new Point2 [n];
// 寻找最近点对
c l o s e ( X , Y, Z , 0 , n - 1 , a , b , d ) ;
// 删除数组并返回
delete [] Y;
delete [] Z;
return true;
}
程序1 4 - 11 计算最近点对
void close(Point1 X[], Point2 Y[], Point2 Z[], int l, int r, Point1& a, Point1& b, float& d)
{//X[l:r] 按x坐标排序
//Y[l:r] 按y坐标排序
if (r-l == 1) {// 两个点
a = X[l];
b = X[r];
d = dist(X[l], X[r]);
r e t u r n ; }
if (r-l == 2) {// 三个点
// 计算所有点对之间的距离
float d1 = dist(X[l], X[l+1]);
float d2 = dist(X[l+1], X[r]);
float d3 = dist(X[l], X[r]);
// 寻找最近点对
if (d1 <= d2 && d1 <= d3) {
a = X[l];
b = X[l+1];
d = d1;
r e t u r n ; }
if (d2 <= d3) {a = X[l+1];
b = X[r];
d = d2;}
else {a = X[l];
b = X[r];
d = d3;}
r e t u r n ; }
/ /多于三个点,划分为两部分
int m = (l+r)/2; // X[l:m] 在A中,余下的在B中
// 在Z[l:m] 和Z [ m + 1 : r ]中创建按y排序的表
int f = l, // Z[l:m]的游标
g = m+1; // Z[m+1:r]的游标
for (int i = l; i <= r; i++)
if (Y[i].p > m) Z[g++] = Y[i];
else Z[f++] = Y[i];
// 对以上两个部分进行求解
c l o s e ( X , Z , Y, l , m , a , b , d ) ;
float dr;
Point1 ar, br;
c l o s e ( X , Z , Y, m + 1 , r, a r, b r, d r ) ;
// (a,b) 是两者中较近的点对
if (dr < d) {a = ar;
b = br;
d = dr;}
M e r g e ( Z , Y,l,m,r);// 重构Y
/ /距离小于d的点放入Z
int k = l; // Z的游标
for (i = l; i <= r; i++)
if (fabs(Y[m].x - Y[i].x) < d) Z[k++] = Y[i];
// 通过检查Z [ l : k - 1 ]中的所有点对,寻找较近的点对
for (i = l; i < k; i++){
for (int j = i+1; j < k && Z[j].y - Z[i].y < d;
j + + ) {
float dp = dist(Z[i], Z[j]);
if (dp < d) {// 较近的点对
d = dp;
a = X[Z[i].p];
b = X[Z[j].p];}
}
}
}
函数c l o s e(见程序1 4 - 11)用来确定X[1:r] 中的最近点对。假定这些点按x 坐标排序。在Y [ 1 : r ]中对这些点按y 坐标排序。Z[ 1 : r ]用来存放中间结果。找到最近点对以后,将在a, b中返回最近点对,在d 中返回距离,数组Y被恢复为输入状态。函数并未修改数组X。
首先考察“小问题”,即少于四个点的点集。因为分割过程不会产生少于两点的数组,因此只需要处理两点和三点的情形。对于这两种情形,可以尝试所有的可能性。当点数超过三个时,通过计算m = ( 1 + r ) / 2把点集分为两组A和B,X [ 1 : m ]属于A,X [ m + 1 : r ]属于B。通过从左至右扫描Y中的点以及确定哪些点属于A,哪些点属于B,可以创建分别与A组和B组对应的,按y 坐标排序的Z [ 1 : m ]和Z [ m + 1 : r ]。此时Y和Z的角色互相交换,依次执行两个递归调用来获取A和B中的最近点对。在两次递归调用返回后,必须保证Z不发生改变,但对Y则无此要求。不过,仅Y [ l : r ]可能会发生改变。通过合并操作(见程序1 4 - 5)可以以Z [ 1 : r ]重构Y [ 1 : r ]。
为实现图1 4 - 1 6的策略,首先扫描Y [ 1 : r ],并收集距分割线小于的点,将这些点存放在Z [ 1 : k - 1 ]中。可按如下两种方式来把RA中点p 与p 的比较区内的所有点进行配对:1) 与RB 中y 坐标≥p.y 的点配对;2) 与y 坐标≤p.y 的点配对。这可以通过将每个点Z [ i ](1≤i < k,不管该点是在RA
还是在RB中)与Z[j] 配对来实现,其中i<j 且Z [ j ] . y - Z [ i ] . y< 。对每一个Z [ i ],在2 × 区域内所检查的点如图1 4 - 1 7所示。由于在每个2 × 子区域内的点至少相距。因此每一个子区域中的点数不会超过四个,所以与Z [ i ]配对的点Z [ j ]最多有七个。
2. 复杂性分析
令t (n) 代表处理n 个点时,函数close 所需要的时间。当n<4时,t (n) 等于某个常数d。当n≥4时,需花费(n) 时间来完成以下工作:将点集划分为两个部分,两次递归调用后重构Y,淘汰距分割线很远的点,寻找更好的第三类点对。两次递归调用需分别耗时t (「n /2ù」和t (?n /2?).
这个递归式与归并排序的递归式完全一样,其结果为t (n) = (nl o gn)。另外,函数c l o s e s t还需耗时(nl o gn)来完成如下额外工作:对X进行排序,创建Y和Z,对Y进行排序。因此分而治之最近点对求解算法的时间复杂性为(nl o gn)。