1. 程式人生 > >分治_求解平面最近點對的O(nlg(n))演算法_演算法分析

分治_求解平面最近點對的O(nlg(n))演算法_演算法分析

          本文主要講述經典的分治法求解平面最近點對的演算法, 為方便敘述, 先給出演算法的C+++實現示例, 然後給出演算法的正確性證明, 並分析其時間複雜度.

//分治法求解平面最近點對
#include <iostream>
#include <cstdio>
#include <vector>
#include <algorithm>
#include <cmath>
using namespace std;
const double NIL = 1e15;
//如果a的縱座標(a.second)小於b的縱座標(b.second)則返回true, 否則返回false 
bool cmpByY(const pair<double, double> a, const pair<double, double> b){
	return a.second < b.second;
}
//返回點a, b之間的距離 
double getDis(const pair<double, double> &a, const pair<double, double> &b){
	return sqrt((a.first - b.first) * (a.first - b.first) 
				+ (a.second - b.second) * (a.second - b.second));
}
//返回點集pointsSet[l...r]中相距最近的兩點之間的距離, 並將pointsSet[l...r]按
//照縱座標(second)非遞減排序, 要求pointsSet[l...r]已經按照橫座標非遞減排序
//, 0 =< l <= r < pointsSet.size();
double getMinDisHelp(vector<pair<double, double> > &pointsSet, int l, int r){
	if(l == r) return NIL;
	int mid = (l + r) >> 1;
	double  midX = pointsSet[mid].first, 
			 ans = min(getMinDisHelp(pointsSet, l, mid), getMinDisHelp(pointsSet, mid + 1, r));
	inplace_merge(&pointsSet[l], &pointsSet[mid + 1], &pointsSet[r + 1], cmpByY);
	vector<pair<double, double> > midPointsSet;
	for(int i = l; i <= r; ++i) 
		if(abs(pointsSet[i].first - midX) <= ans) midPointsSet.push_back(pointsSet[i]); 
	for(int i = 0; i < midPointsSet.size(); ++i)
		for(int j = i + 1; j < midPointsSet.size() && j <= i + 7; ++j)
			ans = min(ans, getDis(midPointsSet[i], midPointsSet[j]));
	return ans;	
}
//返回點集pointsSet中距離最近的兩個點之間的距離 
double getMinDis(const vector<pair<double, double> > &pointsSet){
	vector<pair<double, double> > vectmp = pointsSet;
	return getMinDisHelp(vectmp, 0, vectmp.size() - 1);
}
//測試資料部分 
int main(){
	vector<pair<double, double> > p;
	p.push_back(make_pair(-9, 1)), p.push_back(make_pair(-7, -1)), p.push_back(make_pair(-3, 0))
	, p.push_back(make_pair(-3, 2)), p.push_back(make_pair(-3, -3)), p.push_back(make_pair(1, -1))
	, p.push_back(make_pair(5, -1));
	cout << "PointsSet: (-9, 1), (-7, -1), (-3, 0), (-3, 2), (-3, -3), (1, -1), (5, -1) min distance: " 
		 << getMinDis(p) << endl; 
	return 0;
}

 程式分析: 

對於上述程式, 函式getMinDis(const vector<pair<double, double> > &pointsSet)返回點集pointsSet中距離最近的兩個點之間的距離, 但是我們關注的是其呼叫的輔助函式getMinDisHelp(vector<pair<double, double> > &pointsSet, int l, int r), 該輔助函式返回點集pointsSet[l...r]中相距最近的兩點之間的距離, 並將pointsSet[l...r]按照縱座標(second)非遞減排序, (要求: pointsSet[l...r]已經按照橫座標非遞減排序, 0 =< l <= r < pointsSet.size()).

    定義集合P_{i} = { x | x是一次對於getMinDisHelp(vector<pair<double, double> > &pointsSet, int l, int r)的呼叫, 且滿足r- l + 1 = i, 引數pointsSet[l...r]已經按照橫座標非遞減排序, 0 =< l <= r < pointsSet.size() }

    分析getMinDisHelp的函式體部分, 易知對於集P_{1}, P_{2}中的所有元素均能夠返回點集pointsSet[l...r]中相距最近的兩點之間的距離, 並將pointsSet[l...r]按照縱座標(second)非遞減排序, 為了遞推的考察P_{k}, (k >= 3)

, 需要先解釋一下上述程式碼的第37行為何對區域性計數變數j限制為j <= i + 7, 很容易知道程式的第25行執行完畢後, ans中的值為點對中的點均在直線x = midX同側(可以位於直線x = midX上)的點對的最小距離, 顯而易見, 我們只需要將分佈於直線兩側(具體是一個在pointsSet[l...mid]中, 另一個在pointsSet[mid + 1...r]中)的所有點對距離中的最小值和ans取較小者即為最終答案(集pointsSet[l...r]中相距最近的兩點之間的距離), 顯然如果前者更小, 那麼這樣的點對一定位於直線x = misX - ans和直線x = midX + ans之間(不包括這兩條直線), 對應程式中的minPointsSet. 如下圖所示, 考慮座標系中所有被直線x = midX分成兩個相等的d\times d正方形的2d\times d的矩形區域, 將每個這樣的矩形區域進一步細分為8個相同的邊長為\frac{d}{2}的正方形, 每個正方形對角線長為\frac{\sqrt{2}}{2}d < d, 也就是說位於同一個邊長為\frac{d}{2}的正方形中的兩個點的最大間距小於d, 那麼有結論: 在矩形ACDF中含有的點數不超過8個.下面是此結論的證明:

                                                                 

    很明顯在AKJI, IJLF, OCQN, NQDP這4個正方形(包括邊界)中, 每個正方形內至多含有1個點, 假設在矩形ACDF中含有至少9個點, 那麼在矩形KOPL中至少含有5個點, 不失一般性, 假設有3個點屬於pointsSet[mid + 1...r], 則矩形BOPE中含有至少三個點, 從而正方形BONM或MNPE中至少含有2兩個點, 這兩個點的距離小於d, 與ans為d矛盾, 因此假設不成立, 矩形ACDF中至多含有8個點, 結論得證.

    另一方面, 如果我們將midPoinsSet中的所有點按照縱座標非遞減的方式排序, 那麼位於上述同一個d\times 2d矩形中的點對a, b(不失一般性, 假設a在排序後的序列中位於b之前)之間的所有點也必定位於該矩形區域, 據此可知b一定在a之後的7的點之內, 因此程式第31行對j施加限制j <= i + 7, 並不會導致錯過對最近點對的考察, 在結合第26行的歸併過程便可的結論, 對於任意集P_{i}中的元素均能夠返回點集pointsSet[l...r]中相距最近的兩點之間的距離, 並將pointsSet[l...r]按照縱座標(second)非遞減排序, 從而程式正確性得證

     接下來分析程式的時間複雜度, 顯然我們唯一需要分析的是getMinDisHelp(vector<pair<double, double> > &pointsSet, int l, int r)的時間複雜度, 設n = r - l + 1, 時間複雜度為T(n), 則第15行的兩次遞迴呼叫用時2\times T(\frac{n}{2}), 第26行的歸併和第30行的遍歷均用時O(n), 故可得遞迴式: T(n) = 2\times T(\frac{n}{2}) + O(n), 解此遞迴式得:T(n) = O(nlg(n)), 從而本程式的時間複雜度為O(nlg(n)).