kd树解平面最近点对

 1 #include <cstdio>
 2 #include <cstring>
 3 #include <algorithm>
 4 #include <queue>
 5 #include <cstdio>
 6 #include <cmath>
 7 #define MAX (101000)
 8 #define pow(x) ((x)*(x))
 9 
10 using namespace std;
11 
12 int n, idx;
13 
14 struct Point{
15     double x[2];
16     bool operator < (const Point &u) const{
17         return x[idx] < u.x[idx];
18     }
19 }po[MAX];
20 
21 struct Tree {
22     Point p[MAX<<2];
23     int son[MAX<<2];
24     bool f[MAX<<2];
25     int ps[MAX<<2];
26     void build ( int l , int r , int u = 1, int dep = 0)
27     {
28         if(l > r) return;
29         son[u] = r-l;
30         son[u<<1] = son[u<<1|1] = -1;
31         idx = dep%2;
32         int mid = (l+r)>>1;
33         nth_element(po+l, po+mid, po+r+1);
34         p[u] = po[mid], ps[u] = mid;
35         build ( l , mid-1 , u<<1 , dep+1 );
36         build ( mid+1 , r , u<<1|1 , dep+1 );
37     }
38 
39     double query(Point a, int id, int u = 1, int dep = 0){
40         if(son[u] == -1) return 1000000000.0;
41 
42         double dis = pow(p[u].x[0]-a.x[0]) + pow(p[u].x[1]-a.x[1]);
43         if(ps[u] == id) dis = 1000000000.0;
44         int dim = dep%2, fg = 0;
45         int x = u<<1, y = u<<1|1;
46         if(a.x[dim] >= p[u].x[dim]) swap(x, y);
47         double tm1 = 1000000000.0, tm2 = 1000000000.0;
48         if(~son[x])tm1 = query(a, id, x, dep+1);
49         if(pow(a.x[dim] - p[u].x[dim]) < tm1) fg = 1;
50         if(~son[y] && fg) tm2 = query(a, id, y, dep+1);
51         if(dis > tm1) dis = tm1;
52         if(dis > tm2) dis = tm2;
53         return dis;
54     }
55 }kd;
56 
57 
58 int main(){
59     while(~scanf("%d", &n), n){
60         for(int i = 0; i < n; i++)
61                 scanf("%lf %lf", &po[i].x[0], &po[i].x[1]);
62         memset(kd.f, 0, sizeof(kd.f));
63         memset(kd.ps, -1, sizeof(kd.ps));
64         kd.build(0, n-1);
65         double ans = 1000000000.0;
66         for(int i = 0; i < n; i++){
67             double tm = kd.query(po[i], i);
68             if(tm < ans) ans = tm;
69         }
70         printf("%.2lf\n", sqrt(ans)/2);
71     }
72     return 0;
73 }