LuoguP1742最小圆覆盖 [计算几何]
一道模板题吧,WA了好久
开始以为是 凸包+旋转卡壳+等等 结果才发现有错。
果然模板不可以乱yy
正解:
随机增量法
一看名字就知道,先要把输入的点打乱,使其随机化,作用就是降低复杂度,后面讲。
然后就是从第一个点开始枚举点,如果当前的枚举的点在圆内部,就继续不用管,否者就以该点为圆心半径为0开始枚举前面的点,如果前面的点在当前圆的外面就取点和点的中点为圆心,距离的一半为半径,以这个圆再枚举前面的点,然后如果点在圆外就以三点组成的三角形的外切圆来更新当前的圆。
通过以上操作,因为每次都更新较大的圆,而且答案必定有两个或以上的点在圆上(只有一个点除外),所以保证了正确性,看似其实只有仅仅的复杂度,因为打乱后复杂度平摊下来只有$O(n),证明的具体过程去百度搜吧,应该都比我讲的清楚。
其实还有一道几乎应该就是一样的题,P2533 [AHOI2012]信号塔,这里面的题解有证明。
下面直接上代码
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define db double
const int M=1e5+1;
using namespace std;
const db eps=1e-10;
db dcmp(db x){if(fabs(x)<eps) return 0;return x;}
int n;
struct point{
db x,y;
point(db a=0,db b=0):x(a),y(b){}
void in(){scanf("%lf%lf",&x,&y);}
}pp[M];
point operator +(point a,point b){return point(a.x+b.x,a.y+b.y);}
point operator -(point a,point b){return point(a.x-b.x,a.y-b.y);}
point operator *(point a,db b){return point(a.x*b ,a.y*b );}
point operator /(point a,db b){return point(a.x/b ,a.y/b );}
db cross(point a,point b){return a.x*b.y-a.y*b.x;}
db dot (point a,point b){return a.x*b.x+a.y*b.y;}
double dis(point a,point b)
{return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
point getmid(point a,point b){return point((a.x+b.x)/2,(a.y+b.y)/2);}
point rotate(point a){return point(-a.y,a.x);}
struct circle{
point o;db r;
circle(){}
circle(point a,db r):o(a),r(r){}
};
struct triangle{
point t1,t2,t3;
circle cir2(){
db Bx = t2.x-t1.x, By = t2.y-t1.y;
db Cx = t3.x-t1.x, Cy = t3.y-t1.y;
db D = 2*(Bx*Cy-By*Cx);
db cx = (Cy*(Bx*Bx+By*By) - By*(Cx*Cx+Cy*Cy))/D + t1.x;
db cy = (Bx*(Cx*Cx+Cy*Cy) - Cx*(Bx*Bx+By*By))/D + t1.y;
point p = point(cx, cy);
return circle(p, dis(t1,p));
}//数学方法求外心
triangle(){}
triangle(point a,point b,point c):t1(a),t2(b),t3(c){}
};
circle circlein(point t1,point t2,point t3){
point v1=t2-t1,v2=t1-t3;v1=rotate(v1),v2=rotate(v2);
point p1=getmid(t1,t2),p2=getmid(t1,t3);point u=p1-p2;
db t=cross(v2,u)/cross(v1,v2);
point oo=p1+v1*t;
return circle(oo,dcmp(dis(oo,t1)));
}几何方法求外心,三角形两边中垂线交点
void work(){
circle ans;
ans.o=pp[1];ans.r=0;
for(int i=2;i<=n;i++){
if(dis(pp[i],ans.o)>ans.r+eps){
ans=circle(pp[i],0);
for(int j=1;j<i;j++){
if(dis(pp[j],ans.o)>ans.r+eps){
ans.o=getmid(pp[i],pp[j]);
ans.r=dis(pp[j],pp[i])/2;
for(int k=1;k<j;k++){
if(dis(pp[k],ans.o)>ans.r+eps){
ans=circlein(pp[i],pp[j],pp[k]);
}
}
}
}
}
}
printf("%.10lf\n%.10lf %.10lf\n",ans.r,ans.o.x,ans.o.y);
}
int main()
{
scanf("%d",&n);
for(int i=1;i<=n;i++)pp[i].in();
random_shuffle(pp+1,pp+n+1);//记得打乱
//不打乱的话会被卡复杂度和精度,蒟蒻我被卡了十几次90分。
work();
return 0;
}
Fighting Everyday!