学步园

求点集中两两之间距离之和最小的三个点

在使用增量法过程中，我们需要对于新加入的点是否能够构成新的答案进行高效的判定。假设当前集合的答案为r，对于新加入的点，我们发现如果这个点能够与原集合中另外两点构成新的最近的三个点的话，那么另外的两个点与这个新加入的点的距离一定不超过r/2,证明如下：

假设当前新加入点为C，且C与之前集合中的点A和B构成了新的最近的三个点，如下图所示，不妨设|AC|>=
ans / 2，由于三角形两边之和大于第三边，必有|AB|+ |BC| > |AC|，即
|AB|+ |BC| + |AC| > 2|AC| > ans，也就不可能成为新的答案。

这样我们就可以确定，对于新加入集合的点Pi而言，理论上我们只需要查看以Pi为圆心ans/2为半径的圆内的点即可，但是这样不太容易实现。我们考虑将整个平面划分成变长为ans/2的正方形网格，这样对于Pi而言，我们只需要查找Pi所在网格以及周围8个网格共计9个网格内的点即可。我们可以通过对于网格进行哈希的方法来快速查询这些点。由于随着点的加入，ans越来越小，可以发现任意一个网格中的点总数总是很有限个，这样通过哈希我们就可以再接近常数时间内判断出新加入的点是否能够产生更小的答案。但是这样一来当新加入的点产生了新的答案时，ans发生变化，就要求我们对于更新网格，也就是更新哈希表中的内容，需要O(n)的时间，这样可以得到一个时间复杂度为O(n^2)的增量算法。

但是上述算法仍然不能满足时间要求。我们注意到更新网格不是总会发生，只是最坏情况算法会到O(n^2)的时间复杂度，而实际发生的概率为9k/n，其中k是网格中的点的个数，n为当前集合总的点数，这样我们就可以通过对于点的顺序随机洗牌，然后执行上述增量算法，从而避免最坏情况的发生，也就是使用随机增量法，得到一个O(n)的算法。

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <ctime>
using namespace std;
const int MAXN = 65536;
#define P 971
#define INF 1e20

struct Point {
		double x,y;
}p[MAXN];
struct Hash {
		int id,next;
}hash[MAXN];

int head[MAXN],tot;
const int dx[] = {-1,-1,-1,0,0,0,1,1,1};
const int dy[] = {-1,0,1,-1,0,1,-1,0,1};

double ans;
int cur[MAXN],cs;

inline double dis(const Point &a,const Point &b) {
		return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
inline double cal(const Point &a,const Point &b,const Point &c){
		return dis(a,b) + dis(b,c) + dis(c,a);
}

inline void sherwood(const int n){
		int i,r;
		for(i = 0; i < n; ++i){
				r = rand() % n;
				if (r != i)swap(p[i],p[r]);
		}
}

pair <int,int> getGrid(const Point &p) {
		return make_pair((int)(p.x/(ans/2)+1),(int)(p.y/(ans/2)+1));
}

inline int getKey(const pair<int,int> &grid) {
		return ((grid.first*P) + grid.second)&MAXN;
}

inline void init(){
		memset(head,-1,sizeof(head));
		tot = 0;
}
inline void insert(const int id){
		int i,key = getKey(getGrid(p[id]));
		hash[tot].id = id;
		hash[tot].next = head[key];
		head[key] = tot ++;
}

inline void getPoint(const Point &cp) {
		pair <int,int> grid = getGrid(cp);
		int i,j,key;
		cs = 0;
		for(i = 0; i < 9; ++i){
				key = getKey(make_pair(grid.first+dx[i],grid.second+dy[i]));
				for(j = head[key];~j; j = hash[j].next) {
						if(dis(cp,p[hash[j].id]) < ans/2)cur[cs++] = hash[j].id;
				}
		}
}

inline void rebuild(const int n){
 		int i;
		init();
		for(i = 0; i < n; ++i){
				insert(i);
		}
}
inline void solve(const int n){
		int i,j,k;
		double now;
		sherwood(n);
		ans = cal(p[0],p[1],p[2]);
		rebuild(3);
		for(i = 3; i < n; ++i){
				getPoint(p[i]);
				now = INF;
				for(j = 0; j < cs; ++j){
						if(cur[j] == i)continue;
						for(k = j + 1; k < cs; ++k){
								if(cur[k] == i || cur[k] == cur[j])continue;
								now = min(now,cal(p[i],p[cur[j]],p[cur[k]]));
						}
				}
				if(now < ans){
						ans = now;
						rebuild(i+1);
				}else{
						insert(i);
				}
		}
}


int main(){
		int t,n,i;
		srand((unsigned int)time(NULL));
		scanf("%d",&t);
		while (t --){
				scanf("%d",&n);
				for(i = 0; i < n; ++i){
						scanf("%lf%lf",&p[i].x,&p[i].y);
				}
				solve(n);
				printf("%.3lf\n",ans);
		}
		return 0;
}