LOJ 6409. 「ICPC World Finals 2018」熊猫保护区

官方题解，内附我不会的(O(nlogn))做法
http://www.csc.kth.se/~austrin/icpc/finals2018solutions.pdf
这题题意很简单．

给你一个可以不凸的多边形，要求找一个点距离所有多边形的顶点最小距离最大．

首先这个东西很像是半平面交，因为半平面交可以做

给你一个可以不凸的多边形，要求找一个点距离所有多边形的边最小距离最大．

然而直接上半平面交是搞不出来的．
观察一下性质，显然最优取值点在某一点对的垂直平分线上．

然后大力
https://en.wikipedia.org/wiki/Voronoi_diagram
一发．
发现自己不会写找出维诺图的(O(nlogn))算法．
这时候就有两种思考：
１．乱搞，卡精度
２．观察到数据范围很小，直接暴力半平面交找出维诺图的顶点．（然而我并没有想到这个）

考虑乱搞，怎么乱搞呢，首先要利用结论，垂直平分线上的结论．
那么我随机生成一个点，让他随机沿垂直平分线走一个随机的距离．实践证明，这样做很不理想．
考虑抱精度，那么设一个步长，每次乘上一个alpha，搞一搞精度．
还是不能过．

因为点会走到维诺图外，所以只沿维诺图的某条边所在直线走，不一定可以走到维诺图的最优定点上．
考虑随机一个角度，多随机几次，在那个角度上走步长．

然后发现被近乎一条直线的多边形卡爆了，那么考虑再加一种策略，向一个随机顶点走．
然后还是被卡爆了．

考虑借鉴遗传算法的思想，在下一次的方向中借鉴上一次的较优解．
然后就勉强可以水过．

#include <bits/stdc++.h>

using namespace std;
template<class T>
bool Chkmin(T &x,T y){
	return x>y?(x=y,1):0;
}
template<class T>
bool Chkmax(T &x,T y){
	return x<y?(x=y,1):0;
}
template<class T>
T sqr(T x){
	return x*x;
}
typedef double ld;
typedef long long i64;
const ld INF=1e18;
const ld EPS=1e-8;
struct point{
	ld x,y;
	point operator +(const point &_) const{
		return {x+_.x,y+_.y};
	}
	point operator -(const point &_) const{
		return {x-_.x,y-_.y};
	}
	point operator /(const ld &_){
		return {x/_,y/_};
	}
	point operator *(const ld &_){
		return {x*_,y*_};
	}
	ld cross(const point &_) const{
		return x*_.y-y*_.x;
	}
	ld dot(const point &_) const{
		return x*_.x+y*_.y;
	}
	point rotate() const{
		return {-y,x};
	}
	ld dis(){
		return sqrt(sqr(x)+sqr(y));
	}
};
ostream& operator <<(ostream & os,const point &_){
	return os<<"{"<<_.x<<","<<_.y<<"}";
}
struct segment{
	point s,t;
	bool strict_cross(const segment &_) const{
		return (t-s).cross(_.s-s)*(t-s).cross(_.t-s)<0&&(_.t-_.s).cross(s-_.s)*(_.t-_.s).cross(t-_.s)<0;
	}
	bool operator &(const segment &_) const{
		//cerr<<"&"<<s<<" "<<t<<",,"<<_.s<<" "<<_.t<<endl;
		return strict_cross(_);
	}
};
struct polygon{
	vector<point> data;
	size_t next(size_t val) const{
		return val==size()-1?0:val+1;
	}
	size_t size() const{
		return data.size();
	}
	segment edge(int x) const{
		return segment({data[x],data[next(x)]});
	}
	point& operator [](const size_t &x){
		return data[x];
	}
};
bool in(const point &test,const polygon &con){
	static const point add({42348.23712341,34435.61431236});
	int counter=0;
	for (int i=0; i<con.size(); ++i)
		if (segment({test,test+add})&con.edge(i)) ++counter;
	return counter&1;
}
point neareast(const point &test,polygon &con){
	ld nowdis=(test-con[0]).dis();
	point rvalue=con[0];
	for (size_t i=1; i<con.size(); ++i)
		if (Chkmin(nowdis,(test-con[i]).dis())) rvalue=con[i];
	return rvalue;
}
const int C=10000;
const int T=1000;
point rec[T];
bool b[T];
ld Rand(){
	return (((i64)rand()<<31)+rand());
}
int mxx,mix,mxy,miy;
int RandIntX(){
	return rand()%(mxx-mix+1)+mix;
}
int RandIntY(){
	return rand()%(mxy-miy+1)+miy;
}
polygon con;
point trans(ld x,point y){
	return y*(x/y.dis());
}
point pp(){
	int x,y;
	do{
		x=rand()%con.size();
		y=rand()%con.size();
	}while (x==y);
	return (con[x]-con[y]).rotate();
}
point conclude(int TT){
	point x({RandIntX(),RandIntY()});
	//bug here
	bool fi=1;
	for (ld step=C/2; step>EPS; step*=0.8){
		//cerr<<x<<" "<<step<<endl;
		rec[1]=x+trans(step,x-rec[0]);
		rec[0]=x;
		if (fi)
			for (int i=2,t; i<TT; ++i){
				if ((t=rand())&1){
					rec[i]=x+trans(step,pp());//optimize
				}
				else if ((t>>1)&1) rec[i]=x+trans(step,(con[rand()%con.size()]-x));
				else{
					ld alpha=Rand();
					rec[i]=x+point({step*cos(alpha),step*sin(alpha)});
				}
			}
		else
			for (int i=2,t; i<TT; ++i)
				if ((t=rand())&1){
					ld alpha=Rand();
					rec[i]=x+point({step*cos(alpha),step*sin(alpha)});
				}
				else if ((t>>1)&1) rec[i]=x+trans(step,(con[rand()%con.size()]-x));
				else{
					rec[i]=x+trans(step,pp());//optimize
				}
		int counter=0;
		for (int i=0; i<TT; ++i)
			if (in(rec[i],con)){//bug
				//cerr<<"inint"<<i<<" "<<rec[i]<<endl;
				if (fi){
					//cerr<<x<<endl;
					fi=0;
					step=C/2;
				}
				b[i]=1;
				++counter;
			}
			else b[i]=0;
		//getchar();
		//exit(0);
		if (counter){
			ld nowdis=-1;
			for (int i=0; i<TT; ++i)
				if (b[i]&&Chkmax(nowdis,(rec[i]-neareast(rec[i],con)).dis())){
					//cerr<<"UPD"<<nowdis<<endl;
					x=rec[i];
				}
		}
		else{
			ld nowdis=INF;
			for (int i=0; i<TT; ++i)
				if (Chkmin(nowdis,(rec[i]-neareast(rec[i],con)).dis())){
					x=rec[i];
				}
		}
	}
	return x;
}
int main(){
	int n;
	cin>>n;
	con.data.resize(n);
	mxx=-10000;
	mix=10000;
	mxy=-10000;
	miy=10000;
	for (int i=0; i<n; ++i){
		int x,y;
		cin>>x>>y;
		con[i]=point({x,y});
		mix=min(mix,x);
		mxx=max(mxx,x);
		miy=min(miy,y);
		mxy=max(mxy,y);
	}
	int TIME=clock();
	ld ans=0;
	while (clock()-TIME<CLOCKS_PER_SEC){
		int p=rand()%15+5;
		//cerr<<"P"<<p<<endl;
		point tmp=conclude(p);
		ans=max(ans,(tmp-neareast(tmp,con)).dis());
	}
	while (clock()-TIME<CLOCKS_PER_SEC*1.5){
		int p=rand()%200;
		if (p>100) p=rand()%200;
		if (p>100) p=rand()%200;
		if (p>100) p=rand()%200;
		if (p>100) p=rand()%200;
		p=max(p,30);
		point tmp=conclude(p);
		ans=max(ans,(tmp-neareast(tmp,con)).dis());
		//cerr<<"amns"<<ans<<endl;
	}
	cout<<fixed<<setprecision(10)<<ans;
}

考虑正经算法．
半平面交怎么求维诺图呢？
把某一个顶点作为线段端点的垂直平分线搞出来，然后暴力半平面交，最后的点集就是维诺图上在该点附近的顶点．
结论:维诺图的顶点数和边数都为(O(n))级别．
注意维诺图的边和多边形的交点也有可能是答案．
复杂度(O(n^2logn))

#pragma GCC optimize("Ofast")
#pragma GCC optimize("inline")
#pragma GCC optmize(3)
#include <bits/stdc++.h>
using namespace std;
template<class T>
bool Chkmin(T &x,T y){
	return x>y?(x=y,1):0;
}
template<class T>
bool Chkmax(T &x,T y){
	return x<y?(x=y,1):0;
}
template<class T>
T sqr(T x){
	return x*x;
}
typedef double ld;
typedef long long i64;
const ld INF=1e18;
const ld EPS=1e-9;
struct point{
	ld x,y;
	point operator +(const point &_) const{
		return {x+_.x,y+_.y};
	}
	point operator -(const point &_) const{
		return {x-_.x,y-_.y};
	}
	point operator /(const ld &_) const{
		return {x/_,y/_};
	}
	point operator *(const ld &_) const{
		return {x*_,y*_};
	}
	ld cross(const point &_) const{
		return x*_.y-y*_.x;
	}
	ld dot(const point &_) const{
		return x*_.x+y*_.y;
	}
	bool operator ==(const point &_) const{
		return x==_.x&&y==_.y;
	}
	bool operator !=(const point &_) const{
		return x!=_.x||y!=_.y;
	}
	point rotate() const{
		return {-y,x};
	}
	ld dis() const{
		return sqrt(sqr(x)+sqr(y));
	}
	ld sqrdis() const{
		return x*x+y*y;
	}
	point makelong(ld len=3e4) const{
		ld t=len/(dis());
		return (*this)*t;
	}
};
ld atan2(const point &x){
	return atan2(x.y,x.x);
}
ostream& operator <<(ostream & os,const point &_){
	return os<<"{"<<_.x<<","<<_.y<<"}";
}
struct segment{
	point s,t;
	bool strict_cross(const segment &_) const{
		return (t-s).cross(_.s-s)*(t-s).cross(_.t-s)<0&&(_.t-_.s).cross(s-_.s)*(_.t-_.s).cross(t-_.s)<0;
	}
	bool operator &(const segment &_) const{
		//cerr<<"&"<<s<<" "<<t<<",,"<<_.s<<" "<<_.t<<endl;
		return strict_cross(_);
	}
	bool onright(const point &_) const{
		return (_-s).cross(t-s)>0;
	}
	point generator(const ld &_) const{
		return s+(t-s)*_;
	}
	point intersect(const segment &_) const{
		ld t1=(t-s).cross(_.t-_.s);
		ld t2=(_.t-_.s).cross(s-_.s);
		return generator(t2/t1);
	}
	bool on(const point &_) const{
		//cerr<<(_-s).dot(_-t)<<" "<<fabs((_-s).cross(_-t))<<endl;
		return (_-s).dot(_-t)<0&&fabs((_-s).cross(_-t))<EPS;
	}
	bool line_on(const point &_) const{
		//cerr<<(_-s).dot(_-t)<<" "<<fabs((_-s).cross(_-t))<<endl;
		return (_-s).dot(_-t)<0;
	}
};
ld atan2(const segment &x){
	return atan2(x.t-x.s);
}
ld ans;
struct polygon:vector<point>{
	size_t next(size_t val) const{
		return val==size()-1?0:val+1;
	}
	segment edge(int x) const{
		return segment({(*this)[x],(*this)[next(x)]});
	}
	ld neareast(const point &p) const{
		//cerr<<"ne"<<p<<endl;
		ld ret=INF;
		ld c=ans*ans;
		for (auto i:*this)
			if ((ret=min(ret,(i-p).sqrdis()))<=c) break;
		//if (sqrt(ret)>30000) cerr<<"SUCC"<<p<<" "<<sqrt(ret)<<endl;
		return sqrt(ret);
	}
	bool in(const point &p) const{
		for (size_t i=0; i<this->size(); ++i){
			auto j=edge(i);
			if (j.on(p)) return 1;
		}
		static const point v={43333.3435499546546L,74354.4534534512434L};
		int ret=0;
		for (size_t i=0; i<this->size(); ++i){
			auto j=edge(i);
			ret+=segment{p,p+v}&j;
		}
		return ret&1;
	}
};
typedef vector<segment> vs;
polygon con;
struct Segment:segment{
	ld tt;
	Segment(){
	}
	Segment(const segment &x):segment(x){
		tt=atan2(x);
	}
};
polygon halfplaneintersect(vs v){
	vector<Segment> vp(v.begin(),v.end());
	stable_sort(vp.begin(),vp.end(),[](const Segment &x,const Segment &y)->bool{
			if (x.tt+EPS<y.tt) return 1;
			if (y.tt+EPS<x.tt) return 0;
			return x.onright(y.s);
		});
	//for (auto i:v) cerr<<i.s<<" "<<i.t<<endl;
	//exit(0);
	static const int N=2005;
	static Segment q[N];
	static point p[N];
	int l=0,r=-1;
	for (int i=0; i<vp.size(); ++i){
		//cerr<<"Line: "<<"	"<<v[i].s<<" "<<v[i].t<<" "<<atan2(v[i])<<endl;
		while (l<r&&vp[i].onright(p[r])) --r;
		while (l<r&&vp[i].onright(p[l+1])) ++l;
		if (l<=r&&vp[i].tt<q[r].tt+EPS) continue;
		q[++r]=vp[i];
		if (l<r){
			//cerr<<"intersect 	 	 	"<<l<<" "<<r<<" "<<p[r]<<endl;
			p[r]=q[r-1].intersect(q[r]);
		}
		//cerr<<"this"<<l<<" "<<r<<endl;
	}
	while (l<r&&q[l].onright(p[r])) --r;
	while (l<r&&q[r].onright(p[l+1])) ++l;
	//cerr<<"atlast 	"<<l<<" "<<r<<endl;
	p[l]=q[l].intersect(q[r]);
	polygon ret;
	ret.assign(p+l,p+r+1);
	for (int i=l; i<=r; ++i){
		segment t1=q[i];
		for (int j=0; j<con.size(); ++j){
			segment t2=con.edge(j);
			point tmp=t1.intersect(t2);
			//cerr<<"tmp"<<tmp<<" "<<t1.s<<" "<<t1.t<<" "<<t1.on(tmp)<<endl;
				//cerr<<"init"<<endl;
			if (t2.line_on(tmp))
				ans=max(ans,con.neareast(tmp));
		}
	}
	//for (auto i:ret) cerr<<i<<"|";
	//cerr<<endl;
	return ret;
}
segment perpendicular_bisector(point x,point y){
	point a=(x+y)/2,b=(y-x).rotate().makelong();
	return {a-b,a+b};
}
void solve(polygon result){
	//cerr<<"RE"<<result.size()<<endl;
	for (auto i:result){
		//cerr<<i<<endl;
		if (con.in(i)){
			//cerr<<"I"<<i<<endl;
			ans=max(ans,con.neareast(i));
		}
	}
}
void test(){
	segment a2({{0,0},{1,1}});
	segment a1({{1,0},{-1,0}});
	cerr<<a1.intersect(a2)<<endl;
}
int main(){
	//test();
	int n;
	cin>>n;
	con.resize(n);
	for (int i=0; i<n; ++i){
		int x,y;
		cin>>x>>y;
		con[i]={x,y};
	}
	for (auto i:con){
		//cerr<<"BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB"<<i<<endl;
		vs tmp;
		for (auto j:con)
			if (i!=j){
				tmp.push_back(perpendicular_bisector(i,j));
				//cerr<<tmp.back().s<<" "<<tmp.back().t<<" "<<j<<endl;
			}
		solve(
		halfplaneintersect(tmp)
				);
		//return 0;
	}
	cout<<fixed<<setprecision(10)<<ans;
}