题目链接:
http://acm.hdu.edu.cn/showproblem.php?pid=2966
突然想起来做这个题的原因是早上QQ上一个学弟问了我这么一个问题，平面上有1,000个同半径的圆和20,000个点，问每个圆中包含多少个点。第一个想法是对于每个点，二分一个值k表示有k个圆可以包含这个点，那么这k个圆的圆心一定是距离这个点最近的k个。于是圆心建KD Tree求解。暂时还没想到什么其他高效的做法。
因为KD Tree原来只是听说过从来没接触过，于是去看了下UESTC_木子日匀大神的博文(http://www.mzry1992.com/blog/miao/kd%E6%A0%91.html)和KD Tree原始的Paper Range Searching using KdTre
其实KD Tree的思路非常简单，基本看code很清楚就明白了。这里稍微说一下复杂度分析的过程。
对于一个2维的KD Tree，假设query region中包含$ k $个点，再考虑query region 4条边所在的直线其中的一条$\hat{l}$。假设$\mathbf l(root)$表示划分原始区域的直线，那么$\hat{l}$与$\mathbf l(root)$的左半部分或者右半部分相交。假设与左半部分相交，那么左半部分的两个儿子节点，一定都与$\hat{l}$相交。如果假设我们所求的与$\hat{l}$相交的区域个数为$T(n)$的话，则有
\[T(n) = O(1) \qquad(n = 1) \\ = 2 + 2O(n/4) \qquad(n > 1)\]
根据主定理可得复杂度为$T(n) = O(n^{log_4^2} + k)$ 即 $T(n) = O(\sqrt{n} + k)$
对于查询最近点的复杂度，因为首先对于任意输入点先找到了一个离该点较近的点，然后实质上相当于查询区域选择成为以查询点为圆心当前最优解为半径的圆内的query region。不过这里的$k$到底应该是多少分析不太清楚了..>_<
PS: 为何我用&对齐公式在blog里没有效果，用latex编译通过运行正常...求解
附:
贴个几乎照抄的代码
/*
* Author: fatboy_cw
* Created Time: 2012/12/25 15:01:18
* File Name: hdu2966.cpp
*/
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <cstdlib>
#include <algorithm>
#include <vector>
using namespace std;
#define out(v) cerr << #v << ": " << (v) << endl
#define SZ(v) ((int)(v).size())
const int maxint = -1u>>1;
template <class T> bool get_max(T& a, const T &b) {return b > a? a = b, 1: 0;}
template <class T> bool get_min(T& a, const T &b) {return b < a? a = b, 1: 0;}
typedef long long lint;
const int maxn = 100000 + 5;
struct point {
int x, y;
point() {}
point(int _x, int _y) : x(_x), y(_y) {}
};
bool split[maxn];
point vec[maxn], origin[maxn];
int n, t;
void update(lint &u, lint v) {
if(u == -1 || u > v) {
u = v;
}
}
lint getDist(const point &p1, const point &p2) {
return (lint)(p1.x - p2.x) * (p1.x - p2.x) + (lint)(p1.y - p2.y) * (p1.y - p2.y);
}
bool cmpX(const point &p1, const point &p2) {
return p1.x < p2.x;
}
bool cmpY(const point &p1, const point &p2) {
return p1.y < p2.y;
}
void build(int l, int r, bool spt) {
if(l > r) return;
int mid = (l + r) >> 1;
split[mid] = spt;
if(!spt) nth_element(vec + l, vec + mid, vec + r + 1, cmpX);
else nth_element(vec + l, vec + mid, vec + r + 1, cmpY);
build(l, mid - 1, !spt);
build(mid + 1, r, !spt);
}
void find(int l, int r, const point &p, lint &res) {
if(l > r) return;
int mid = (l + r) >> 1;
lint t_res = getDist(p, vec[mid]);
if(t_res) update(res, t_res);
lint dis = split[mid] ? (p.y - vec[mid].y) : (p.x - vec[mid].x);
int l1 = l, r1 = mid - 1, l2 = mid + 1, r2 = r;
if(dis > 0) swap(l1, l2), swap(r1, r2);
find(l1, r1, p, res);
if(res == -1 || dis * dis < res) {
find(l2, r2, p, res);
}
}
int main() {
scanf("%d", &t);
while(t--) {
scanf("%d", &n);
for(int i = 0; i < n; i++) {
int x, y;
scanf("%d%d", &x, &y);
vec[i] = origin[i] = point(x, y);
}
build(0, n - 1, 0);
for(int i = 0; i < n; i++) {
lint res = -1;
find(0, n - 1, origin[i], res);
printf("%I64d\n", res);
}
}
return 0;
}
