Peienwu's Blog

Hello world

0%

最近點對：詳解4種不同複雜度之算法

發表於 分類於 ， Valine：
文章字數： 8.4k 所需閱讀時間 ≈ 8 分鐘

題目敘述
給定二維平面上 $n$ 個點，每一點都有座標 $(x_i,y_i)$ ，求出最近的點對之歐幾里德距離為多少？
$dis(p_i,p_j) = \sqrt{(x_1-x_2)^2+(y_1-y_2)^2}$

平面最近點對有好多種實作方式，從最差的暴力枚舉、稍微優化的掃描線演算法、到分治與隨機，有4種不同的時間複雜度。利用TIOJ 1500這一題最近點對的裸題，來實測各種不同複雜度下所需要的執行時間。

暴力枚舉

時間複雜度：$O(N^2)$

Submission
時間：TLE,10440

暴力$O(n^2)$將所有點進行枚舉，因為值域是 $n≤50000$ ，平方枚舉會有TLE的問題。

程式碼

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
#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

這一種作法是改善過後的暴力枚舉，利用計算幾何中掃描線的概念，先將所有點依照x座標進行排序（y座標隨意）。接著想像一條從左往右掃的掃描線，對於每一條掃描線看右邊的點，如果當前最近點對距離為 $d$，因此只要遇上x座標差距大於 $d$ 的點時，即可繼續下一輪的枚舉。

加上排序的關係，其時間複雜度至少為 $O(n\log n)$，但這種掃描線的方式無法有效過濾所有點都在相同的x座標上的情況，因此最差的時間複雜度會退化成 $O(n^2)$ ，不過聽說平均的狀況下是很快的！

上圖為掃描線執行最近點對的一個示意圖，黑線為掃描線，$d$ 為掃描線左邊所有點的最近點對距離，我們只要每一輪枚舉這個點與右邊x座標差在 $d$ 以內的所有點，即可進行下一輪的更新！

程式碼

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
#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

原本以為上面的掃描線就是他的極限了，沒想到上面的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中。

程式碼

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
#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

分治做最近點對的基本想法，先將所有點依照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)$。（不過實際上應該會比這個快，因為並不是要對所有點都進行排序）。

如果要做得更快，可以在y座標排序的地方稍微動動手腳。既然每一層都要對y座標進行排序，排序好的東西再排序一次其實沒有什麼意義，因此就可以用）合併排序（merge sort）的方式，將所有已經排序好的兩個左右序列進行$O(n)$的合併（可以用std::merge()完成），如此一來，就不須要每一層花到 $O(n\log n)$ 的ㄕˊ間進行排序，使總複雜度降低為 $O(n\log n)$！

程式碼

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
#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

用隨機算法做最近點對的期望複雜度是 $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)$。

程式碼

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
#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;
    }
}

相關題目