[计算几何]平面最近点对

【算法分析】:

暴力枚举O(N^2)

分治解决O(NlogN)

传送门:http://blogold.chinaunix.net/u2/63316/showart_2334130.html

没有多想……按照算法的原理自己写了一个……效率还是蛮高的!10W的数据秒杀,暴力则跑了2min多种才出解。

正确性无从而知,因此拿暴力对拍了一下,还好自测数据都是对的。

【Code】:

const
	maxn=100001;
type
	point=record
		x,y:extended;
	end;
var
	q,q1,q2:array[0..maxn]of longint;
	p,pp:array[0..maxn]of point;
	n,i:longint;
function min(x,y:extended):extended;
	begin
		if x<y then exit(x); exit(y);
	end;
procedure swap(var a,b:point);
	var t:point;
	begin
		t:=a; a:=b; b:=t;
	end;
procedure qsort(l,r:longint);
	var i,j:longint; x:extended;
	begin
		i:=l; j:=r; x:=p[(l+r)>>1].x;
		repeat
			while p[i].x<x do inc(i);
			while p[j].x>x do dec(j);
			if i<=j then begin
				swap(p[i],p[j]);
				inc(i); dec(j);
			end;
		until i>j;
		if l<j then qsort(l,j);
		if i<r then qsort(i,r);
	end;
function dist(p0,p1:point):extended;
	begin
		exit(sqrt(sqr(p0.x-p1.x)+sqr(p0.y-p1.y)));
	end;
function solve(l,r:longint):extended;
	var tl,tr,tx,mx:extended; m,i,j,k:longint;
	begin
		if l=r then exit(1e99);                                //区间大小为1
		if l+1=r then begin
			if p[l].y>p[r].y then swap(p[l],p[r]);
			exit(dist(p[l],p[r]));
		end;                                                   //区间大小为2
		m:=(l+r)>>1; mx:=(p[m].x+p[m+1].x)/2;
		tl:=solve(l,m);                                        //左区间
		tr:=solve(m+1,r);                                      //右区间
		tx:=min(tl,tr);
		q1[0]:=0; q2[0]:=0; q[0]:=0;
		for i:=l to m do
			if p[i].x>=mx-tx then begin
				inc(q1[0]); q1[q1[0]]:=i;
			end;                                               //筛出左边可能的点
		for i:=m+1 to r do
			if p[i].y<=mx+tx then begin
				inc(q2[0]); q2[q2[0]]:=i;
			end;                                               //筛出右边肯能的点
		j:=1; k:=1;
		for i:=1 to q1[0]+q2[0] do begin
			if (j>q1[0]) then q[i]:=q2[k] else
			if (k>q2[0]) then q[i]:=q1[j] else
			if (p[q1[j]].y<p[q2[k]].y) then q[i]:=q1[j] else q[i]:=q2[k];
			if (q[i]=q1[j]) and (j<=q1[0]) then inc(j) else inc(k);
		end;                                                   //合并到一个线性表中
		for i:=1 to q1[0]+q2[0] do
			for j:=1 to 6 do if (i+j<=q1[0]+q2[0]) then
				tx:=min(tx,dist(p[q[i]],p[q[i+j]]));           //枚举周围6个的可能的解
		j:=l; k:=m+1;
		for i:=l to r do begin
			if (j>m) then pp[i]:=p[k] else
			if (k>r) then pp[i]:=p[j] else
			if (p[j].y<p[k].y) then pp[i]:=p[j] else pp[i]:=p[k];
			if (pp[i].x=p[j].x) and (pp[i].y=p[j].y)  and (j<=m) then inc(j) else inc(k);
		end;
		for i:=l to r do p[i]:=pp[i];                          //将区间按照y坐标排序
		exit(tx);
	end;
begin
	readln(n);
	for i:=1 to n do readln(p[i].x,p[i].y);
	qsort(1,n);                                                //将区间按照x坐标排序
	writeln(solve(1,n):0:4);
end.

posted on 2011-03-17 11:22  Skywalker_Q  阅读(344)  评论(0)    收藏  举报

导航