平面最近點對:比較4種不同複雜度之算法


原文網址:https://peienwu.com/pair/

問題定義:給定二維平面上 $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座標在範圍內的點進行枚舉更新答案。

  1. 將點輸入並且排序,X座標為主,Y座標為輔。
  2. 使用set,並以Y座標為排序基準(pair的首項),以儲存第 $i$ 點的左方、水平距離小於等於d的點。
  3. 右掃描線依序窮舉各點作為右端點。
     (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有做好的話,期望可以在線性時間解決這個問題。基本的想法如下:

  1. 將最近點對距離設為d,初始為第一、二個點之間的距離
  2. 將每一個點的座標塞入以 $\frac{d}{2}$ 為邊長的網格中
  3. 將點加入網格中,查看要加入的網格是否已經有點在其中
  4. 一個網格不可容納兩個點,否則必須更新最近點對的距離
  5. 在更新最近點對距離之後,將前面的點的網格座標以新的$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;
    }
}

相關題目

#程式設計 #最近點對 #closest pair problem #分治 #掃描線 #隨機






你可能感興趣的文章

第八週作業訂正、Twitch API

第八週作業訂正、Twitch API

用 Python 自學資料科學與機器學習入門實戰:Numpy 基礎入門

用 Python 自學資料科學與機器學習入門實戰:Numpy 基礎入門

氣泡排序(Bubble Sort)、插入排序(Insertion Sort)、選擇排序(Selection Sort)

氣泡排序(Bubble Sort)、插入排序(Insertion Sort)、選擇排序(Selection Sort)






留言討論