LuoguP1742最小圆覆盖 [计算几何]

一道模板题吧,WA了好久

开始以为是 凸包+旋转卡壳+等等 结果才发现有错。
果然模板不可以乱yy

正解:

随机增量法

一看名字就知道,先要把输入的点打乱,使其随机化,作用就是降低复杂度,后面讲。

然后就是从第一个点开始枚举点i,如果当前的枚举的点在圆内部,就继续不用管,否者就以该点为圆心半径为0开始枚举i前面的点j,如果前面的点在当前圆的外面就取点i和点j的中点为圆心,距离的一半为半径,以这个圆再枚举j前面的点k,然后如果点k在圆外就以i,j,k三点组成的三角形的外切圆来更新当前的圆。

通过以上操作,因为每次都更新较大的圆,而且答案必定有两个或以上的点在圆上(只有一个点除外),所以保证了正确性,看似O(n3)其实只有仅仅O(n)的复杂度,因为打乱后复杂度平摊下来只有$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!

posted @ 2018-06-11 19:24  VictoryCzt  阅读(135)  评论(0编辑  收藏  举报