容易想到二分实数,关键是如何判定一个半径是否合法。
这个判定问题其实可以转化为判定能否将多边形的边分成两部分,每一部分都可以用多边形内的一点作半径为的圆与其内部所有边所在直线相交。给定了某一部分的边的话,将该部分的边向内缩半径后可以得到若干半平面,多边形的边也对应了若干半平面,显然合法的点需要在这些半平面内,问题即为判定这些半平面与多边形对应的半平面的交是否为空。
注意到若某个集合的半平面与多边形对应的半平面的交不为空,一定包含某两条直线(包括向内缩的和原多边形的)的交点。一个暴力的算法是枚举所有直线的交点(在多边形内的),再枚举每条边向内缩后对应的半平面看是否在内部,这样得到一个可以覆盖的边的集合,最后看能不能用两个这样的点覆盖所有边,时间复杂度为,不能通过。
要得到复杂度更优秀的算法,我们需要证明一个结论,将边按逆时针顺序标号后,最优解中分成的两部分都是连续的区间(可以循环,即这样的,下同)。
CF上原题解的证明过程有问题,这里感谢毛啸同学(mathew99)给出的证明。
设我们最后选择的点为和,一般地,只考虑和不重合的情况。恰当地旋转坐标系和选择原点,可以使和落在轴上且关于原点对称(令在原点左侧)。考虑凸多边形每条边所在直线与轴的交点(若与轴平行的话与和距离相等,可以忽略),那么它不可能落在和之间,否则与凸多边形的定义矛盾,并且显然会划分成连续的两部分,一部分交点落在左侧,一部分落在右侧(事实上每一部分的交点坐标是单峰的)。注意到若交点落在左侧,会距离更近;若落在右侧,会距离更近,因此距离更近的和距离更近的边都形成连续的区间(相等的划分到一侧即可),最优解中显然也会这么划分。于是证明完毕。
有了这个结论后,我们的想法是算出每个区间能否用一个点作圆覆盖。如果直接枚举每个区间计算半平面交判定,复杂度并没有优化。不过我们可以注意到若区间的子区间不能用一个点作圆覆盖,那也不能用一个点作圆覆盖。于是我们可以考虑双指针,枚举左端点,那么最大可行的单调不降,因此总共只用判定个区间。算出了每个最大可行的后就容易判定是否可行了。
实现的时候我们需要预先对所有的半平面排序,再利用增量法求出半平面交。时间复杂度。
#include <bits/stdc++.h>
#define eps 1e-10
using namespace std;
typedef double db;
inline int sgn(db x) {
return (x>-eps)?((x>eps)?1:0):-1;
}
struct Point {
db x,y;
Point() {}
Point(db a,db b):x(a),y(b) {}
db ang() {return atan2(y,x);}
db len() {return sqrt(x*x+y*y);}
Point nor() {
db l=len();
return Point(-y/l,x/l);
}
Point operator + (Point b) {return Point(x+b.x,y+b.y);}
Point operator - (Point b) {return Point(x-b.x,y-b.y);}
Point operator * (db b) {return Point(x*b,y*b);}
Point operator / (db b) {return Point(x/b,y/b);}
};
inline db cross(Point x,Point y) {
return x.x*y.y-x.y*y.x;
}
struct Line {
Point x,y;
db ang;
Line() {}
Line(Point a,Point b):x(a),y(b) {ang=(y-x).ang();}
bool operator < (const Line & b) const {
if (sgn(ang-b.ang)) return ang<b.ang;
Point t1=y,t2=b.y;
return sgn(cross(t1-x,t2-x))<0;
}
};
inline Point inter(Line x,Line y) {
db s1=cross(y.x-x.x,x.y-x.x);
db s2=cross(x.y-x.x,y.y-x.x);
return (y.x*s2+y.y*s1)/(s1+s2);
}
Line p[605];
bool in[605];
Point curp;
bool half_plane(int n) {
static Line q1[605];
static Point q2[605];
int l=1,r=0;
for(int i=1;i<=n;i++)
if (in[i]&&(l>r||sgn(p[i].ang-q1[r].ang))) {
Line t=p[i];
while (l<r&&sgn(cross(t.y-t.x,q2[r]-t.x))<0) r--;
while (l<r&&sgn(cross(t.y-t.x,q2[l+1]-t.x))<0) l++;
q1[++r]=t;
if (l<r) q2[r]=inter(q1[r-1],t);
}
while (l<r&&sgn(cross(q1[l].y-q1[l].x,q2[r]-q1[l].x))<0) r--;
if (r-l<2) return 0;
curp=q2[r];
return 1;
}
Line a[305],b[305],q[605];
int id1[305],id2[305],rpos[305];
Point bel[305];
bool cmp(int x,int y) {
return q[x]<q[y];
}
bool check_half(int n,int l,int r) {
memset(in,0,sizeof(in));
for(int i=1;i<=n;i++) in[id1[i]]=1;
for(int i=l;i<=r;i++) in[id2[(i-1)%n+1]]=1;
return half_plane(2*n);
}
Point ans1,ans2;
bool check(int n,db r) {
static int cur[605];
for(int i=1;i<=n;i++) {
Point t=(a[i].y-a[i].x).nor();
b[i]=Line(a[i].y+t*r,a[i].x+t*r);
}
for(int i=1;i<=n;i++) {
q[i]=a[i];q[n+i]=b[i];
}
for(int i=1;i<=2*n;i++) cur[i]=i;
sort(cur+1,cur+2*n+1,cmp);
for(int i=1;i<=2*n;i++) {
if (cur[i]<=n) id1[cur[i]]=i;
else id2[cur[i]-n]=i;
p[i]=q[cur[i]];
}
int nr=0;
for(int i=1;i<=n;i++) {
while (nr<2*n&&check_half(n,i,nr+1)) nr++;
check_half(n,i,nr);
rpos[i]=nr;bel[i]=curp;
}
if (rpos[1]>=n) {
ans1=ans2=bel[1];
return 1;
}
for(int i=1;i<n;i++)
for(int j=i+1;j<=n;j++)
if (rpos[i]>=j-1&&rpos[j]>=n+i-1) {
ans1=bel[i];ans2=bel[j];
return 1;
}
return 0;
}
Point fir[305];
int main() {
int n;
scanf("%d",&n);
for(int i=1;i<=n;i++) {
int x,y;
scanf("%d%d",&x,&y);
fir[i]=Point(x,y);
}
for(int i=1;i<=n;i++) a[i]=Line(fir[i],fir[i%n+1]);
db l=0,r=1e5;
while (r-l>1e-7) {
db mid=(l+r)*.5;
if (check(n,mid)) r=mid; else l=mid;
}
check(n,l+5e-8);
printf("%.10f\n",l);
printf("%.10f %.10f\n",ans1.x,ans1.y);
printf("%.10f %.10f\n",ans2.x,ans2.y);
return 0;
}