【BZOJ4219】跑得比谁都快（Voronoi图）（Delaunay三角剖分）

传送门

---

D剖板子题，参考了这篇文章的讲解：here，以及darklove的实现：here

确实是细节一大堆，很多地方能不能取等号要好好斟酌一下。太多了不想讲

目前找得到的C++代码只有darklove的是严格$O(n\log n)$，其他都是$O(n^2)$没有被卡。。。

---

代码：

#include<bits/stdc++.h>
#define ll long long
#define re register
#define db double
#define cs const

namespace IO{
	inline char gc(){
		static cs int Rlen=1<<22|1;static char buf[Rlen],*nl,*nr;
		return (nl==nr)&&(nr=(nl=buf)+fread(buf,1,Rlen,stdin),nl==nr)?EOF:*nl++;
	}template<typename T>T get(){
		char c;bool f=false;while(!isdigit(c=gc()))f=c=='-';T num=c^48;
		while(isdigit(c=gc()))num=((num+(num<<2))<<1)+(c^48);return f?-num:num;
	}inline int gi(){return get<int>();}
	template<typename T>T get_float(){
		char c;bool f=false;while(!isdigit(c=gc()))f=c=='-';T x=c^48;
		while(isdigit(c=gc()))x=x*10+(c^48);if(c!='.')return f?-x:x;
		T y=1.0;while(isdigit(c=gc()))x+=(y/=10)*(c^48);return f?-x:x;
	}inline db gd(){return get_float<db>();}
}using namespace IO;

using std::cerr;
using std::cout;

cs int N=1e5+7;

cs db eps=1e-10;

struct Pnt{
	db x,y;Pnt(){}Pnt(db _x,db _y):x(_x),y(_y){}
	friend Pnt operator+(cs Pnt &a,cs Pnt &b){return Pnt(a.x+b.x,a.y+b.y);}
	friend Pnt operator-(cs Pnt &a,cs Pnt &b){return Pnt(a.x-b.x,a.y-b.y);}
	friend db operator*(cs Pnt &a,cs Pnt &b){return a.x*b.y-b.x*a.y;}
	friend db dot(cs Pnt &a,cs Pnt &b){return a.x*b.x+a.y*b.y;}
	db len()cs{return sqrt(x*x+y*y);}
};

struct pt{
	db x,y,z;pt(){}pt(db _x,db _y,db _z):x(_x),y(_y),z(_z){}
	friend pt operator+(cs pt &a,cs pt &b){return pt(a.x+b.x,a.y+b.y,a.z+b.z);}
	friend pt operator-(cs pt &a,cs pt &b){return pt(a.x-b.x,a.y-b.y,a.z-b.z);}
	friend pt operator*(cs pt &a,cs pt &b){return pt(a.y*b.z-b.y*a.z,a.z*b.x-b.z*a.x,a.x*b.y-a.y*b.x);}
	friend db dot(cs pt &a,cs pt &b){return a.x*b.x+a.y*b.y+a.z*b.z;}
	pt operator-()cs{return pt(-x,-y,-z);}
	db len()cs{return sqrt(x*x+y*y+z*z);}
};

inline db crs(cs Pnt &a,cs Pnt &b,cs Pnt &c){
	return (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x);
}
inline pt crs(cs pt &a,cs pt &b,cs pt &c){
	return (b-a)*(c-a);
}
inline pt trans(cs Pnt &a){
	return pt(a.x,a.y,a.x*a.x+a.y*a.y);
}
inline bool incir(cs Pnt &D,cs Pnt &A,cs Pnt &B,cs Pnt &C){
	pt a=trans(A),b=trans(B),c=trans(C),d=trans(D);
	pt dir=crs(a,b,c);if(dir.z>0)dir=-dir;return dot(d-a,dir)>eps;
}

cs int M=N<<3|7;

int n;db X,Y;

Pnt p[N];
int bin[M],tot;
int pr[M],nx[M];
int to[M],rev[M];

inline int ins_pr(int u,int v){
	int e=bin[++tot];
	nx[e]=u,pr[e]=pr[u];
	nx[pr[e]]=e,pr[u]=e;
	to[e]=v;return e;
}
inline int ins_nx(int u,int v){
	int e=bin[++tot];
	pr[e]=u,nx[e]=nx[u];
	pr[nx[e]]=e,nx[u]=e;
	to[e]=v;return e;
}
inline void del(int e){
	pr[nx[e]]=pr[e];
	nx[pr[e]]=nx[e];
	bin[tot--]=e;
}

inline db crs(int a,int b,int c){
	return crs(p[a],p[b],p[c]);
}
inline bool incir(int d,int a,int b,int c){
	return incir(p[d],p[a],p[b],p[c]);
}

void delaunay(int l,int r){
	if(l==r)return ;
	if(l+1==r){
		int el=ins_nx(l,r);
		int er=ins_pr(r,l);
		rev[el]=er,rev[er]=el;
		return ;
	}int mid=(l+r)>>1,tp=0;static int st[N];
	delaunay(l,mid),delaunay(mid+1,r);
	for(int re i=l;i<=r;++i){
		while(tp>1&&crs(st[tp-1],st[tp],i)<-eps)--tp;
		st[++tp]=i;
	}int nl=0,nr=0;
	for(int re i=1;i<tp;++i)
		if(st[i]<=mid&&st[i+1]>mid)
			{nl=st[i],nr=st[i+1];break;}
	int sl,sr;
	for(sl=pr[nl];sl!=nl;sl=pr[sl])
		if(crs(nl,nr,to[sl])>-eps)break;
	for(sr=nx[nr];sr!=nr;sr=nx[sr])
		if(crs(nr,nl,to[sr])<eps)break;
	int el=ins_nx(sl,nr),er=ins_pr(sr,nl);
	rev[el]=er,rev[er]=el;
	while(true){
		bool cl=false,cr=false;
		for(;sl!=nl;sl=pr[sl]){
			if(crs(nl,nr,to[sl])<eps){
				if(p[to[sl]].x>=p[nl].x)continue;
				break;
			}
			if(pr[sl]!=nl&&incir(to[pr[sl]],nl,nr,to[sl]))
				del(rev[sl]),del(sl);
			else {cl=true;break;}
		}
		for(;sr!=nr;sr=nx[sr]){
			if(crs(nr,nl,to[sr])>-eps){
				if(p[to[sr]].x<p[nr].x)continue;
				break;
			}
			if(nx[sr]!=nr&&incir(to[nx[sr]],nr,nl,to[sr]))
				del(rev[sr]),del(sr);
			else {cr=true;break;}
		}
		if(!cl&&!cr)break;
		if(!cr||(cl&&!incir(to[sr],nl,nr,to[sl]))){
			nl=to[sl];
			for(sl=pr[nl];sl!=nl;sl=pr[sl])
				if(crs(nl,nr,to[sl])>-eps||p[to[sl]].x<p[nl].x)break;
			el=ins_nx(sl,nr),er=ins_pr(sr,nl);
			rev[el]=er,rev[er]=el;
		}else {
			nr=to[sr];
			for(sr=nx[nr];sr!=nr;sr=nx[sr])
				if(crs(nr,nl,to[sr])<eps||p[to[sr]].x>p[nr].x)break;
			el=ins_nx(sl,nr),er=ins_pr(sr,nl);
			rev[el]=er,rev[er]=el;
		}
	}
}

struct edge{int u,v;db w;};
std::vector<edge> E;int fa[N];
inline int gf(int u){while(u^fa[u])u=fa[u]=fa[fa[u]];return u;}

void Main(){
	n=gi(),X=gd(),Y=gd();
	for(int re i=1;i<=n;++i)p[i].x=gd(),p[i].y=gd();
	tot=n;for(int re i=1;i<M;++i)bin[i]=i;
	for(int re i=1;i<=n;++i)pr[i]=nx[i]=i;
	std::sort(p+1,p+n+1,[](cs Pnt &a,cs Pnt &b){
		return a.x<b.x||(a.x==b.x&&a.y<b.y);
	});delaunay(1,n);for(int re i=1;i<=n;++i){
		E.push_back({i,n+1,std::min(p[i].x,Y-p[i].y)});
		E.push_back({i,n+2,std::min(p[i].y,X-p[i].x)});
		for(int re e=nx[i];e!=i;e=nx[e])if(to[e]>i)
			E.push_back({i,to[e],(p[i]-p[to[e]]).len()*0.5});
	}
	std::sort(E.begin(),E.end(),
		[](cs edge &a,cs edge &b){return a.w<b.w;});
	for(int re i=1;i<=n+2;++i)fa[i]=i;
	for(auto &e:E)if(gf(e.u)!=gf(e.v)){
		fa[gf(e.u)]=gf(e.v);
		if(gf(n+1)==gf(n+2)){
			printf("%.6lf",e.w);return ;
		}
	}
}

void file(){
#ifdef zxyoi
	freopen("faster.in","r",stdin);
#endif
}
signed main(){file();Main();return 0;}