closest pairs

(Updated: )

closest pairs

divide and conquer

对于所有点按 坐标排序，然后左右分开分治

假设左边包含 的点，右边包含 的点，那么中间点设为 号点

两边会分别对点按照 坐标排序，合并时要先归并排序两边的所有点

假设两边计算后最小距离为 , 每次合并时，只考虑 坐标距离中间点 以内的点

此时这些点按照 坐标排好序

枚举每个点，我们同样只考虑距离当前点 坐标不超过 的点

这样的点是有限个的，所以这步直接暴力枚举，也可以双指针

怎么证明有限个，感性理解一下，在一个 的长方形内，所有点距离不小于 ， 每落一个点相当于挖去半径为 的圆，那么有限步后整个长方形就空了

也可以理解为用圆覆盖长方形，有限步后就全覆盖了

最后复杂度

code:

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
#include<bits/stdc++.h>
#define int long long 
using namespace std;
const int N=4e5+10,INF=1e15+10;
int T,n;
pair<int, int> p[N], q[N], t[N];
int dis(pair<int,int> l, pair<int,int> r){
    return (l.first-r.first)*(l.first-r.first)+
        (l.second-r.second)*(l.second-r.second);
}
int divide(int l, int r){
    if(r - l <= 1){
        if(l == r) return INF;
        else{
            if(p[l].second>p[r].second) swap(p[l], p[r]);
            return dis(p[l], p[r]);
        }
    }
    int mid=(l+r)>>1;
    pair<int, int> middle = p[mid+1];
    int d=min(divide(l,mid),divide(mid+1,r));
    int tot=0,ll=l,rr=mid+1,cnt=0;
    while(ll<=mid && rr<=r){
        if(p[ll].second<p[rr].second) q[++tot] = p[ll++];
        else q[++tot] = p[rr++];
    }
    while(ll<=mid) q[++tot] = p[ll++];
    while(rr<=r) q[++tot] = p[rr++];
    for(int i=l;i<=r;++i) p[i]=q[i-l+1];
    for(int i=1;i<=tot;++i){
        if((q[i].first-middle.first)*(q[i].first-middle.first)<=d) t[++cnt]=q[i];
    }
    for(int i=1;i<=cnt;++i){
        for(int j=i-1;j>=1 && (t[j].second-t[i].second)*(t[j].second-t[i].second)<=d;--j){
            d=min(d,dis(t[i],t[j]));
        }
    }
    return d;
}
bool cmp(pair<int,int> x, pair<int,int> y){
    return x.first<y.first;
}
void solve() {
    scanf("%lld",&n);
    for(int i=1,x,y;i<=n;++i){
        scanf("%lld%lld",&x,&y);
        p[i]=make_pair(x, y);
    }
    sort(p+1,p+n+1,cmp);
    printf("%lld\n",divide(1,n));
}
signed main() {
    solve();
    return 0;
}

sieve closest-pair algorithm

一个期望 的神奇筛算法，不是乱搞

第一部分

我们考虑当前点的集合

从中随机挑选一个节点 , 计算出 中其他点到他的最小距离

令 , 创建一个边长为 的网格，假设点 的横纵坐标为 , 那么每个点处在 对应的网格中

定义一个点的邻居为自己所在方格加八邻接的八个方格中的点的集合，记为 。

那么我们将 的孤立点全部删去

重复上面步骤直至集合为空，记最后一次的 , 则得到的最小距离为

第二部分

以这个距离为边长做出网格，计算所有点到他邻居的距离，这些距离的最小值即为答案

分析

Lemma 0: 第一部分每次循环会将所有 的点删去，所有 的点留下

一个点距离邻居的最远距离不超过 (跨过两个网格的对角线), 距离非邻居的最近距离不小于 (直跨过一个网格)

  • Fact 1: 如果一个点距离另一个点距离至少为 , 那他们一定互相不是邻居

  • Fact 2: 如果一个点距离另一个点距离至多为 , 那他们一定互相是邻居

可以改为

因此第一部分每次循环会将所有 的点删去，所有 的点留下

可能会被删去

Lemma 1: 第一部分结束时，满足 , 表示全局最优解（所有点最近距离）

首先 是最优解，所以一定满足

假设第一个使得最近点对中的点被删掉的循环为第 个

首先我们有所有 的点都会被删去，所以 单调不增

假设这个最近点对中的点为 ，那么有 , 同时我们有所有 的点都不会被删去，所以 一定满足

所以

Lemma 2: , 即第一部分的开销是 的

在第一部分中，我们定义了一个点的邻居

第一部分的求孤立点只需要查看自己和周围八个网格中的点数量即可，这一步可以通过把网格哈希做到

所以每一个循环花费是 的

由 Fact 1, 2 我们知道第一部分每次循环会将所有 的点删去，所有 的点留下

考虑按照 将点排序

由于我们随机取点，所以期望每次选的点是中间点，那么点右侧的所有点会被删去，而左侧一部分点也可能删去，所以平均每次会有 至少一半 的点被删去

那么假设衰减系数为 , 那么等比求和，, 所以花费是 的

Lemma 3: 第二部分的复杂度是 的

假设网格边长 满足 , 那么我们有：

  • Fact 3: 每个点的邻居是有限个
  • Fact 4: 每个最近点对中的点都存在于另一个点的邻居中

Fact 3 是由于任意点对的距离 , 那么用半径为 的圆覆盖，有限步就会把九个网格全覆盖

Fact 4 是由于 Fact 1, 的点一定是邻居, 而这保证了我们计算所有点的邻居点可以找到答案

所以这部分的复杂度是 的

那么全部算法期望复杂度是 的

code:

用 map 替代 hash:

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
#include<bits/stdc++.h>
#define int long long 
using namespace std;
const int N=2e6+10,INF=1e15+10;
int T,n;
typedef pair<int, int> pi; 

map<pi, int> vis;
map<pi, vector<pi> > f;
vector<pi> p, q, t;
int dis(pi l, pi r){
    return (l.first-r.first)*(l.first-r.first)+
        (l.second-r.second)*(l.second-r.second);
}
bool cmp(pi x, pi y){
    return x.first<y.first;
}
bool check(int x, int y){ // is isolate
    for(int i=x-1;i<=x+1;++i){
        for(int j=y-1;j<=y+1;++j){
            pi now = make_pair(i, j);
            if(vis[now]>(i==x && j==y)) return false;
        }
    }
    return true;
}
int calculate(pi me, int x, int y){
    int d = INF;
    for(int i=x-1;i<=x+1;++i){
        for(int j=y-1;j<=y+1;++j){
            pi now = make_pair(i, j);
            for(auto you : f[now]){
                if(me != you) d = min(d, dis(me, you));
            }
        }
    }
    return d;
}
void work(){
    mt19937 rnd(time(0));
    vector<int> D;
    t = p;
    while(!p.empty()){
        int size = p.size(), id = rnd()%size;
        int d = INF;
        vis.clear();
        for(int i=0;i<size;++i){
            if(i == id) continue;
            d = min(d, dis(p[i], p[id]));
        }
        double b = sqrt(d)/3.0;
        q.clear();
        for(int i=0;i<size;++i){
            int x = (int)(p[i].first/b), y = (int)(p[i].second/b);
            vis[make_pair(x, y)]++;
        }

        for(int i=0;i<size;++i){
            int x = (int)(p[i].first/b), y = (int)(p[i].second/b);
            if(!check(x, y)) q.push_back(p[i]);
        }
        p = q;
        D.push_back(d);
    }
    double d = sqrt(D[D.size()-1]);
    int size = t.size();
    for(int i=0;i<size;++i){
        int x=(int)(t[i].first/d),y=(int)(t[i].second/d);
        f[make_pair(x, y)].push_back(t[i]);
    }
    int ans=INF;
    for(int i=0;i<size;++i){
        int x=(int)(t[i].first/d),y=(int)(t[i].second/d);
        ans=min(ans,calculate(t[i], x, y));
    }
    printf("%lld\n", ans);
}
void solve() {
    scanf("%lld",&n);
    for(int i=1,x,y;i<=n;++i){
        scanf("%lld%lld",&x,&y);
        p.push_back(make_pair(x, y));
    }
    work();
}

signed main() {
    solve();
    return 0;
}
//g++ -std=c++17 P7883_On.cpp -o main.exe

这个代码 P1429 和 P7883 都过不了，会超时

hash版本：

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
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
#include<bits/stdc++.h>
#define int long long 
using namespace std;
const int N=1e7+30,INF=2e18+10,mod=1e7+19,MOD=1e9+7,B=131,b=233;
int T,n;
typedef pair<int, int> pi; 

int vis[N];
vector<pi> f[N];
//map<pi, vector<pi> > f;
vector<pi> p, q, t;
int dis(pi l, pi r){
    return (l.first-r.first)*(l.first-r.first)+
        (l.second-r.second)*(l.second-r.second);
}
bool cmp(pi x, pi y){
    return x.first<y.first;
}
int hashing(int x, int y){
    int res = 0; 
    x = (x * MOD) % mod + mod;
    y = (y * MOD) % mod + mod;
    res = ((x * B + y * b) % MOD + MOD) % MOD;
    int ans = res % mod;
    return ans;
}
bool check(int x, int y){ // is isolate
    for(int i=x-1;i<=x+1;++i){
        for(int j=y-1;j<=y+1;++j){
            int now = hashing(i, j);
            if(vis[now]>(i==x && j==y)) return false;
        }
    }
    return true;
}

int calculate(pi me, int x, int y){
    int d = INF;
    for(int i=x-1;i<=x+1;++i){
        for(int j=y-1;j<=y+1;++j){
            int now = hashing(i, j);
            for(auto you : f[now]){
                d = min(d, dis(me, you));
            }
        }
    }
    return d;
}
void work(){
    mt19937 rnd(time(0));
    vector<int> D;
    t = p;
    while(!p.empty()){
        int size = p.size(), id = rnd()%size;
        int d = INF;
        for(int i=0;i<size;++i){
            if(i == id) continue;
            d = min(d, dis(p[i], p[id]));
        }
        double b = sqrt(d)/3.0;
        q.clear();
        for(int i=0;i<size;++i){
            int x = (int)(p[i].first/b), y = (int)(p[i].second/b);
            vis[hashing(x, y)]++;
        }
        for(int i=0;i<size;++i){
            int x = (int)(p[i].first/b), y = (int)(p[i].second/b);
            if(!check(x, y)) q.push_back(p[i]);
        }
        for(int i=0;i<N;++i) vis[i]=0; // this is faster lol
        // for(int i=0;i<size;++i){
        //     int x = (int)(p[i].first/b), y = (int)(p[i].second/b);
        //     vis[hashing(x, y)]--;
        // }
        p = q;
        D.push_back(d);
    }
    double d = sqrt(D[D.size()-1]);
    if(d == 0){
        puts("0.0000");
        return;
    }
    int size = t.size();
    int ans = INF;
    for(int i=0;i<size;++i){
        int x=(int)(t[i].first/d),y=(int)(t[i].second/d);
        ans = min(ans,calculate(t[i], x, y));
        f[hashing(x, y)].push_back(t[i]);
    }
    printf("%lld\n", ans);
    //printf("%.4lf\n", sqrt(ans));
}
void solve() {
    scanf("%lld",&n);
    for(int i=1,x,y;i<=n;++i){
        scanf("%lld%lld",&x,&y);
        p.push_back(make_pair(x, y));
    }
    work();
}

signed main() {
    solve(); 
    return 0;
}
/*
g++ -std=c++17 P7883_On_hashing.cpp -o main.exe
main.exe
*/

这里这个 hash 函数并不好，会有很多冲突，结果就是 P1429 可以较大概率过 (冲突了可能会死循环)， P7883 稳定过不了，主要是因为计算邻居和 hash 函数常数大，以及 hash 函数冲突导致第一部分每次循环不能平均去掉一半以上的点 (有其他 hash 值相同但实际不相邻的点使得这个点不是孤立的)