官方题解，内附我不会的(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;
}