問題定義:給定二維平面上 $n$ 個點,每一點都有座標 $(x_i,y_i)$ ,求出最近的點對之歐幾里德距離為多少?
定義 $d(p_i,p_j) = \sqrt{(x_1-x_2)^2+(y_1-y_2)^2}$
最近點對有好多種實作方式,從最差的暴力枚舉、稍微優化的掃描線演算法、到分治與隨機,有4種不同的時間複雜度。利用TIOJ 1500這一題最近點對的裸題,來實測各種不同複雜度下所需要的執行時間。
暴力枚舉
時間複雜度:$O(N^2)$
Submission
時間:TLE,10440(ms)
暴力$O(n^2)$將所有點進行枚舉,因為值域是 $n≤50000$ ,平方枚舉會有TLE的問題。
程式碼
#include <bits/stdc++.h>
#define Orz ios::sync_with_stdio(0),cin.tie(0)
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define pii pair<int,int>
#define pdd pair<double,double>
#define int long long
#define ll long long
#define ld long double
#define N 50005
#define all(x) x.begin(),x.end()
#define INF 5e18
#define eps 1e-9
#define x first
#define y second
using namespace std;
int n;
pii p[N];
ld dis(pii a, pii b){
ld x = a.x-b.x, y = a.y-b.y;
return sqrt(x*x + y*y);
}
signed main(){
Orz;
cout<<fixed<<setprecision(6);
while(cin>>n){
rep(i,0,n-1)cin>>p[i].x>>p[i].y;
ld d = INF;
rep(i,0,n-1){
rep(j,i+1,n-1){
d = min(d, dis(p[i],p[j]));
}
}
cout<<d<<endl;
}
}
掃描線算法
時間複雜度:Worst Case $O(N^2)$
Submission
時間:AC,1668(ms)
這一種作法是改善過後的暴力枚舉,利用計算幾何中掃描線的概念,先將所有點依照x座標進行排序(y座標隨意)。接著想像一條從左往右掃的掃描線,對於每一條掃描線看右邊的點,如果當前最近點對距離為 $d$,因此只要遇上x座標差距大於 $d$ 的點時,即可繼續下一輪的枚舉。
加上排序的關係,其時間複雜度至少為 $O(n\log n)$,但這種掃描線的方式無法有效過濾所有點都在相同的x座標上的情況,因此最差的時間複雜度會退化成 $O(n^2)$ ,不過聽說平均的狀況下是很快的!
上圖為掃描線執行最近點對的一個示意圖,黑線為掃描線,$d$ 為掃描線左邊所有點的最近點對距離,我們只要每一輪枚舉這個點與右邊x座標差在 $d$ 以內的所有點,即可進行下一輪的更新!
程式碼
#include <bits/stdc++.h>
#define Orz ios::sync_with_stdio(0),cin.tie(0)
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define pii pair<int,int>
#define pdd pair<double,double>
#define int long long
#define ll long long
#define ld long double
#define N 50005
#define all(x) x.begin(),x.end()
#define INF 5e18
#define eps 1e-9
#define x first
#define y second
using namespace std;
int n;
pii p[N];
ld dis(pii a, pii b){
ld x = a.x-b.x, y = a.y-b.y;
return sqrt(x*x + y*y);
}
signed main(){
Orz;
cout<<fixed<<setprecision(6);
while(cin>>n){
rep(i,0,n-1)cin>>p[i].x>>p[i].y;
sort(p,p+n);
ld d = INF;
rep(i,0,n-1){
rep(j,i+1,n-1){
if(p[j].x > p[i].x + d)break;
d = min(d, dis(p[i],p[j]));
}
}
cout<<d<<endl;
}
}
掃描線算法(優化後)
時間複雜度:$O(N\log N)$
Submission
時間:AC,148(ms)
原本以為上面的掃描線就是他的極限了,沒想到上面的worst case還可以透過set優化成 $O(n\log n)$!簡單來說,方法一樣是想像一條掃描線由左而右,一樣照上面的想法,把x座標差大於d的點排除,之後利用set二分搜找出y座標在範圍內的點進行枚舉更新答案。
- 將點輸入並且排序,X座標為主,Y座標為輔。
- 使用set,並以Y座標為排序基準(pair的首項),以儲存第 $i$ 點的左方、水平距離小於等於d的點。
- 右掃描線依序窮舉各點作為右端點。
(1) Erase與右端點水平距離大於d的點們(左掃描線右移)
(2) 用二分搜找出與第 $i$ 點垂直距離小於d的點,並嘗試更新
(3) 將第 $i$ 點加入set中。
程式碼
#include <bits/stdc++.h>
#define int long long
#define ld long double
#define N 200005
#define x first
#define y second
#define pii pair<int,int>
#define IOS ios::sync_with_stdio(0),cin.tie(0)
using namespace std;
int n;
vector<pii> p;
set<pii> s;
ld dis(pii a, pii b){
ld x = a.x-b.x, y = a.y-b.y;
return sqrt(x*x + y*y);
}
signed main(){
IOS;
cout<<fixed<<setprecision(6);
while(cin>>n){
p.assign(n,{0,0});
for(int i = 0;i < n;i++)cin>>p[i].x>>p[i].y;
sort(p.begin(),p.end());
s.clear();
s.insert({p[0].y,p[0].x});
int l = 0;ld ans = 5e18;
for(int i = 1;i < n;i++){
int d = ceil(ans);
while(l < i && p[l].x < p[i].x - d){
s.erase({p[l].y,p[l].x});
l++;
}
auto it_l = s.lower_bound({p[i].y - d,0});
auto it_r = s.upper_bound({p[i].y + d,0});
for(auto it = it_l;it != it_r;it++){
ans = min(ans,dis({it->y,it->x},p[i]));
}
s.insert({p[i].y,p[i].x});
}
cout<<ans<<endl;
}
}
分治算法
時間複雜度:$O(N\log N)$
Submission
時間:AC,196(ms)
分治做最近點對的基本想法,先將所有點依照x座標排序,利用遞迴得到分割點左右兩邊所有點的最短距離(兩點並不會跨過中間分隔線),枚舉所有會橫跨兩側且有可能更新最短距離的點對。
從兩半邊的遞迴得到目前的最近點對距離 $d = min(d_l,d_r)$ ,將分隔線附近x座標差距小於$d$的點通通都枚舉一遍。可能會有一個疑問,我們是不是可以縮小枚舉的範圍,否則點的數量可能會太多導致複雜度爆炸?除了x座標可以做點的篩選之外,在枚舉的過程中,我們會利用將所有點對y座標排序,將y座標直線距離大於 $d$ 的情況剔除,所剩下真的需要枚舉點也只會剩下常數個,因此可以放心枚舉。
複雜度分析: 腦海中想像遞迴樹的長相,會發現每一層都需要都需要對y座標進行排序,時間為$O(n\log n)$ ,每一次都將n的值除以2,因此共有$O(\log n)$ 層,總共的時間複雜度為 $O(n\log^2n)$。(不過實際上應該會比這個快,因為並不是要對所有點都進行排序)。
$$T(n) = 2T(\frac{n}{2})+O(n\log n) = O(n\log^2n)$$
如果要做得更快,可以在y座標排序的地方稍微動動手腳。既然每一層都要對y座標進行排序,排序好的東西再排序一次其實沒有什麼意義,因此就可以用)合併排序(merge sort)的方式,將所有已經排序好的兩個左右序列進行$O(n)$的合併(可以用std::merge()完成),如此一來,就不須要每一層花到 $O(n\log n)$ 的ㄕˊ間進行排序,使總複雜度降低為 $O(n\log n)$!
$$T(n) = 2T(\frac{n}{2})+O(n) = O(n\log n)$$
程式碼
#include <bits/stdc++.h>
#define ll long long
#define int long long
#define double long double
#define Orz ios::sync_with_stdio(0),cin.tie(0)
#define N 50002
#define INF1 100000000
#define INF 5e18
#define FOR(i,n) for(int i=0;i<n;i++)
#define rep(i,l,r) for(int i=l;i<=r;i++)
#define pii pair<int,int>
#define x first
#define y second
#define pid pair<int,double>
#define pdi pair<double,int>
#define pdd pair<double,double>
using namespace std;
int n;
vector<pii> p,temp;
void init(){
cout<<fixed<<setprecision(6);
temp.clear();
p.assign(n,{0,0});
}
bool cmp(pii a,pii b){
return a.y < b.y;
}
double dis(pii a,pii b){
double x1 = a.x-b.x,y1 = a.y-b.y;
return sqrt(x1 * x1 + y1 * y1);
}
//區間[l,r]
double solve(int l,int r){
if(l == r)return INF;
int mid = (l+r)/2,mid_pos = p[mid].x;;
double ans = min(solve(l,mid),solve(mid+1,r));
temp.assign((r-l+1),{0,0});
merge(
p.begin() + l, p.begin() + mid + 1,
p.begin() + mid + 1, p.begin() + r + 1,
temp.begin(), cmp
);
rep(i, l, r)p[i] = temp[i-l];
temp.clear();
rep(i, l, r){
if(abs(p[i].x - mid_pos) <= ans){
temp.push_back(p[i]);
}
}
int len = temp.size();
rep(i, 0, len-1){
rep(j, i+1, len-1){
ans = min(ans, dis(temp[i],temp[j]));
if(abs(temp[i].y-temp[j].y) > ans)
break;
}
}
return ans;
}
signed main(){
Orz;
while(cin>>n){
init();
rep(i,0,n-1)cin>>p[i].x>>p[i].y;
sort(p.begin(),p.end());
double ans = solve(0,n-1);
cout<<ans<<endl;
}
}
隨機算法
時間複雜度:期望 $O(N)$
Submission
時間:AC,488(ms)
用隨機算法做最近點對的期望複雜度是 $O(n)$ ,也就是說如果一開始進行的Random_shuffle有做好的話,期望可以在線性時間解決這個問題。基本的想法如下:
- 將最近點對距離設為d,初始為第一、二個點之間的距離
- 將每一個點的座標塞入以 $\frac{d}{2}$ 為邊長的網格中
- 將點加入網格中,查看要加入的網格是否已經有點在其中
- 一個網格不可容納兩個點,否則必須更新最近點對的距離
- 在更新最近點對距離之後,將前面的點的網格座標以新的$d$進行更新
這個算法用到隨機的因子,因此如果在一開始有將所有點進行均勻的打散的話,可以做到期望複雜度 $O(n)$。
複雜度分析:
考慮加入第i+1個點時出現新的最近點對,發生的機率為:在$C_2^{i+1}$個配對中跟i+1個點產生最近點對共有i種可能因此機率為$\frac{2}{i+1}$。
當機率發生的時候,必須將所有的點都刪掉重新來一遍(r變小,重新推入i+1個點),需要付出$O(i+1)$的時間,相乘起來加入每一個點期望的複雜度為$O(1)$,因此總時間複雜度為$O(n)$。
程式碼
#include <bits/stdc++.h>
#define int long long int
#define ld long double
#define ios ios::sync_with_stdio(0),cin.tie(0)
#define N 200005
#define INF 1000000000LL
#define swift 1000000000
using namespace std;
int n;
ld r,d,ans;
int dx[25] = {-2,-2,-2,-2,-2,-1,-1,-1,-1,-1,0,0,0,0,0,1,1,1,1,1,2,2,2,2,2};
int dy[25] = {-2,-1,0,1,2,-2,-1,0,1,2,-2,-1,0,1,2,-2,-1,0,1,2,-2,-1,0,1,2};
unordered_map<int, int> m;
void solve();
inline void init();
void solve();
bool insert(int,int,int);
inline double dis(int,int);
inline int Grid(int);
struct node{
int x,y,ind;
}point[N];
//函式實作
inline void init(){
m.clear();
cout<<fixed<<setprecision(6);
}
inline int Grid(int ind){ //input網格座標
int x = point[ind].x/r;
int y = point[ind].y/r;
return x*INF+y;
}
inline ld dis(node a,node b){
ld x = a.x-b.x,y = a.y-b.y;
return sqrt(x*x+y*y);
}
void solve(){
m.insert(make_pair(Grid(0),0));m.insert(make_pair(Grid(1),1));
for(int ind = 2;ind < n;ind++){
int x = point[ind].x/r,y = point[ind].y/r,better=0;
for(int i=0;i<25;i++){
int nx = x+dx[i],ny = y+dy[i];
auto it = m.find(nx*INF+ny);
if(it!=m.end()){
double distance = dis(point[it->second],point[ind]);
if(distance<d){
better = 1;
ans = dis(point[it->second],point[ind]);
d = distance;
r = d/2;
}
}
}
if(better){
m.clear();
for(int i=0;i<=ind;i++)m.insert(make_pair(Grid(i),i));
}
else{
m.insert(make_pair(Grid(ind), ind));
}
}
}
signed main(){
ios;
while(cin>>n){
init();
for(int i=0;i<n;i++){
int x,y;cin>>x>>y;
x+=swift;y+=swift;
point[i].x = x;point[i].y = y;
}
random_shuffle(point, point+n);
int smalln = sqrt(n);
ans = dis(point[0],point[1]);
d = dis(point[0], point[1]);
for(int i=0;i<=smalln;i++){
for(int j=i+1;j<=smalln;j++){
d = min(d,dis(point[i], point[j]));
ans = min(ans,dis(point[i],point[j]));
}
}
r = d/2;
solve();
cout<<ans<<endl;
}
}