计算几何全家桶

#include<bits/stdc++.h>
/*一:【准备工作】*/
#define db double
#define ll long long
#define Vector Point
using namespace std;
const db eps=1e-8,Pi=acos(-1.0);
inline int dcmp(db a){return a<-eps?-1:(a>eps?1:0);}//处理精度
inline db Abs(db a){return a*dcmp(a);}//取绝对值
struct Point
{
    db x,y;Point(db X=0,db Y=0){x=X,y=Y;}
    void In(){cin>>x>>y;}
    void Out(){cout<<fixed<<setprecision(10)<<x<<" "<<y<<endl;}
};
/*二:【向量】*/
db Dot(Vector a,Vector b){return a.x*b.x+a.y*b.y;}//【点积】
db Cro(Vector a,Vector b){return a.x*b.y-a.y*b.x;}//【叉积】
db Len(Vector a){return sqrt(Dot(a,a));}//【模长】
db Len2(Vector a){return Dot(a,a);} 
db Angle(Vector a,Vector b){return acos(Dot(a,b)/Len(a)/Len(b));}//【两向量夹角】
Vector Normal(Vector a){return Vector(-a.y,a.x);}//【法向量】
Vector operator+(Vector a,Vector b){return Vector(a.x+b.x,a.y+b.y);}
Vector operator-(Vector a,Vector b){return Vector(a.x-b.x,a.y-b.y);}
Vector operator*(Vector a,db b){return Vector(a.x*b,a.y*b);}
bool operator==(Point a,Point b){return !dcmp(a.x-b.x)&&!dcmp(a.y-b.y);}//两点坐标重合则相等

/*三:【点、向量的位置变换】*/

/*1.【点、向量的旋转】*/
Point turn_P(Point a,db theta){//【点A\向量A顺时针旋转theta(弧度)】
    db x=a.x*cos(theta)+a.y*sin(theta);
    db y=-a.x*sin(theta)+a.y*cos(theta);
    return Point(x,y);
}
Point turn_PP(Point a,Point b,db theta){//【将点A绕点B顺时针旋转theta(弧度)】
    db x=(a.x-b.x)*cos(theta)+(a.y-b.y)*sin(theta)+b.x;
    db y=-(a.x-b.x)*sin(theta)+(a.y-b.y)*cos(theta)+b.y;
    return Point(x,y);
}

/*四:【图形与图形之间的关系】*/

/*1.【点与线段】*/
int OnSeg(Point p,Point a,Point b){//【判断点P是否在线段AB上】
    return !dcmp(Cro(p-a,b-a))&&dcmp(Dot(p-a,p-b))<=0;//做法一
//  return !dcmp(Cro(p-a,b-a))&&dcmp(min(a.x,b.x)-p.x)<=0&&dcmp(p.x-max(a.x,b.x))<=0&&dcmp(min(a.y,b.y)-p.y)<=0&&dcmp(p.y-max(a.y,b.y))<=0;//做法二
    //PA,AB共线且P在AB之间(其实也可以用len(p-a)+len(p-b)==len(a-b)判断,但是精度损失较大)
}
db dis_PS(Point p,Point a,Point b){//【点P到线段AB距离】
    if(a==b)return Len(p-a);//AB重合
    Vector x=p-a,y=p-b,z=b-a;
    if(dcmp(Dot(x,z))<0)return Len(x);//P距离A更近
    if(dcmp(Dot(y,z))>0)return Len(y);//P距离B更近
    return Abs(Cro(x,z)/Len(z));//面积除以底边长
}

/*2.【点与直线】*/
int OnLine(Point p,Point a,Point b){//【判断点P是否在直线AB上】
    return !dcmp(Cro(p-a,b-a));//PA,AB共线
}
db dis_PL(Point p,Point a,Point b){//【点P到直线AB距离】
    if(a==b)return Len(p-a);//AB重合
    return Abs(Cro(p-a,b-a)/Len(b-a));//面积除以底边长
}
Point FootPoint(Point p,Point a,Point b){//【点P到直线AB的垂足】
    Vector x=p-a,y=p-b,z=b-a;
    db len1=Dot(x,z)/Len(z),len2=-1.0*Dot(y,z)/Len(z);//分别计算AP,BP在AB,BA上的投影
    return a+z*(len1/(len1+len2));//点A加上向量AF
}
Point Symmetry_PL(Point p,Point a,Point b){//【点P关于直线AB的对称点】
    return p+(FootPoint(p,a,b)-p)*2;//将PF延长一倍即可
}

/*3.【线与线】*/
Point cross_LL(Point a,Point b,Point c,Point d){//【两直线AB,CD的交点】
    Vector x=b-a,y=d-c,z=a-c;
    return a+x*(Cro(y,z)/Cro(x,y));//点A加上向量AF
}
int pancross_LS(Point a,Point b,Point c,Point d){//【判断直线AB与线段CD是否相交】
    return OnSeg(cross_LL(a,b,c,d),c,d);//直线AB与直线CD的交点在线段CD上
}
int pancross_SS(Point a,Point b,Point c,Point d){//【判断两线段AB,CD是否相交】
    if(OnSeg(a,c,d)||OnSeg(b,c,d)||OnSeg(c,a,b)||OnSeg(d,a,b))return 1; 
    db c1=Cro(b-a,c-a),c2=Cro(b-a,d-a);
    db d1=Cro(d-c,a-c),d2=Cro(d-c,b-c);
    return dcmp(c1)*dcmp(c2)<0&&dcmp(d1)*dcmp(d2)<0;//分别在两侧
}
db dis_SS(Point a,Point b,Point c,Point d)//【线段AB,CD的距离】
{
	if(pancross_SS(a,b,c,d))return 0;
	db s1=dis_PS(a,c,d),s2=dis_PS(b,c,d),s3=dis_PS(c,a,b),s4=dis_PS(d,a,b); 
	return min(s1,min(s2,min(s3,s4))); 
}
int Relation_LL(Point a,Point b,Point c,Point d){//【判断直线AB与直线CD的关系】 
    if(!dcmp(Cro(b-a,b-c)))
    {
    	if(!OnLine(a,c,d))return 1;//平行
		return 2;//重合 
	}
	if(!dcmp(Dot(b-a,d-c)))return 3;//垂直 
    return 0; 
}

/*4.【多边形】*/
db PolyLensum(Point *P,int n){//【任意多边形P的周长】
    db S=0;P[n+1]=P[1]; 
    for(int i=1;i<=n;i++)S+=Len(P[i+1]-P[i]);
    return S;
}
db PolyArea(Point *P,int n){//【任意多边形P的面积】
    db S=0;P[n+1]=P[1];
    for(int i=1;i<=n;i++)S+=Cro(P[i],P[i+1]);
    return S/2.0;
}
int PolyShape(Point *P,int n){//【判断任意多边形P的凹凸】
    int op=0;P[0]=P[n];P[n+1]=P[1]; 
	for(int i=1;i<=n;i++)
	{
		int o=dcmp(Cro(P[i+1]-P[i],P[i]-P[i-1]));
		if(!o)continue;
		if(!op)op=o;
		else if(op!=o)return 0;//凹多边形 
	}
	return 1;//凸多边形
}
/*5.【点与多边形】*/
int PIP(Point *P,int n,Point a){//【射线法】判断点A是否在任意多边形Poly以内
    int cnt=0;db tmp;
    for(int i=1;i<=n;++i){
        int j=i<n?i+1:1;
        if(a==P[i]||a==P[j])return 3;//点在定点上 
        if(OnSeg(a,P[i],P[j]))return 2;//点在多边形上
        if(a.y>=min(P[i].y,P[j].y)&&a.y<max(P[i].y,P[j].y))//纵坐标在该线段两端点之间
            tmp=P[i].x+(a.y-P[i].y)/(P[j].y-P[i].y)*(P[j].x-P[i].x),cnt+=dcmp(tmp-a.x)>0;//交点在A右方
    }
    return cnt&1;//穿过奇数次则在多边形以内
}
int judge(Point a,Point L,Point R){//判断AL是否在AR右边
    return dcmp(Cro(L-a,R-a))>0;//必须严格以内
}
int PIP_(Point *P,int n,Point a){//【二分法】判断点A是否在凸多边形Poly以内
    //点按逆时针给出
    if(judge(P[1],a,P[2])||judge(P[1],P[n],a))return 0;//在P[1_2]或P[1_n]外
    if(OnSeg(a,P[1],P[2])||OnSeg(a,P[1],P[n]))return 2;//在P[1_2]或P[1_n]上
    int l=2,r=n-1;
    while(l<r){//二分找到一个位置pos使得P[1]_A在P[1_pos],P[1_(pos+1)]之间
        int mid=(l+r)/2+1;
        if(judge(P[1],P[mid],a))l=mid;
        else r=mid-1;
    }
    if(judge(P[l],a,P[l+1]))return 0;//在P[pos_(pos+1)]外
    if(OnSeg(a,P[l],P[l+1]))return 2;//在P[pos_(pos+1)]上
    return 1;
}
/*6.【线与多边形】*/

/*7.【多边形与多边形】*/
int judge_PP(Point *A,int n,Point *B,int m){//【判断多边形A与多边形B是否相离】
    for(int i1=1;i1<=n;++i1){
        int j1=i1<n?i1+1:1;
        for(int i2=1;i2<=m;++i2){
            int j2=i2<m?i2+1:1;
            if(pancross_SS(A[i1],A[j1],B[i2],B[j2]))return 0;//两线段相交
            if(PIP(B,m,A[i1])||PIP(A,n,B[i2]))return 0;//点包含在内
        }
    }
    return 1;
}
/*六:【凸包】*/
/*1.【求凸包】*/
bool cmp1(Point a,Point b){return a.x!=b.x?a.x<b.x:a.y<b.y;};//按坐标排序
int ConvexHull(Point *P,int n,Point *tb){//【水平序Graham扫描法(Andrew算法)】求凸包
    sort(P+1,P+n+1,cmp1);int t=0;
    for(int i=1;i<=n;i++){//下凸包
        while(t>1&&dcmp(Cro(tb[t]-tb[t-1],P[i]-tb[t-1]))<=0)--t;
        tb[++t]=P[i];
    }
    int St=t;
    for(int i=n-1;i>0;i--){//上凸包
        while(t>St&&dcmp(Cro(tb[t]-tb[t-1],P[i]-tb[t-1]))<=0)--t;
        tb[++t]=P[i];
    }
    return --t;//第一个点算了两次,要减一
}
/*2.【旋转卡壳】*/
int Calipers(Point *p,int n)
{
	//点按逆时针给出 
	db ans=0;p[n+1]=p[1];
	for(int i=1,j=2;i<=n;i++)
	{
		while(dcmp(Cro(p[i+1]-p[i],p[j]-p[i])-Cro(p[i+1]-p[i],p[j+1]-p[i]))<0)j=j<n?j+1:1;
		ans=max(ans,max(Len2(p[j]-p[i]),Len2(p[j]-p[i+1]))); 
	}
	return ans;
}
/*3.【半平面交】*/
struct Line{Point a,b;};
db Atan(Vector a){return atan2(a.y,a.x);}
Point Cross(Line L1,Line L2){return cross_LL(L1.a,L1.b,L2.a,L2.b);}//获取直线L1,L2的交点
int judge(Line L,Point a){return dcmp(Cro(a-L.a,L.b-L.a))>=0;}//判断点a是否在直线L的右边(在线上也要算)
bool cmp2(Line L1,Line L2)
{
	db t1=Atan(L1.b-L1.a),t2=Atan(L2.b-L2.a);
	if(dcmp(t1-t2)!=0)return t1<t2;
	return dcmp(Cro(L2.a-L1.a,L1.b-L1.a))>0;
}
int Halfcut(Line *L,int n,Point *P){//【半平面交】
    sort(L+1,L+n+1,cmp2);int h=1,t=0,m=n;Line Q[n+1];n=0;
    for(int i=1;i<=m;i++)if(i==1||dcmp(Atan(L[i].b-L[i].a)-Atan(L[i-1].b-L[i-1].a)))L[++n]=L[i];
    for(int i=1;i<=n;i++)
	{
        while(h<t&&judge(L[i],Cross(Q[t],Q[t-1])))t--;//当队尾两个直线交点不是在直线L[i]上或者左边时就出队
        while(h<t&&judge(L[i],Cross(Q[h],Q[h+1])))h++;//当队头两个直线交点不是在直线L[i]上或者左边时就出队
        Q[++t]=L[i];
    }
    while(h<t&&judge(Q[h],Cross(Q[t],Q[t-1])))t--;
    while(h<t&&judge(Q[t],Cross(Q[h],Q[h+1])))h++;
    m=0;
    if(h<t)for(int i=h;i<=t;++i)P[++m]=Cross(Q[i],Q[i<t?i+1:h]);
    return m;
}
/*4.【闵可夫斯基和】*/
int Mincowski(Point *P1,int n,Point *P2,int m,Vector *V)
{//【闵可夫斯基和】求两个凸包{P1},{P2}的向量集合{V}={P1+P2}构成的凸包
    Vector V1[n+1],V2[m+1];
    for(int i=1;i<=n;++i)V1[i]=P1[i<n?i+1:1]-P1[i];
    for(int i=1;i<=m;++i)V2[i]=P2[i<m?i+1:1]-P2[i];
    int t=0,i=1,j=1;V[++t]=P1[1]+P2[1];
    while(i<=n&&j<=m)t++,V[t]=V[t-1]+(Cro(V1[i],V2[j])>0?V1[i++]:V2[j++]);
    while(i<=n&&j<=m)t++,V[t]=V[t-1]+(Cro(V1[i],V2[j])>0?V1[i++]:V2[j++]);
    while(i<=n)t++,V[t]=V[t-1]+V1[i++];
    while(j<=m)t++,V[t]=V[t-1]+V2[j++];
    return t;
}
/*七:【圆】*/

/*1.【三点确定一圆】*/
struct Circle{Point O;db r;Circle(Point P={0,0},db R=0){O=P,r=R;}};
Circle GetCircle(Point A,Point B,Point C){//【三点确定一圆】向量垂心法
    Point P1=(A+B)*0.5,P2=(A+C)*0.5;
    Point O=cross_LL(P1,P1+Normal(B-A),P2,P2+Normal(C-A));
    return Circle(O,Len(A-O));
}

/*2.【最小覆盖圆】*/
int PIC(Circle C,Point a){return dcmp(Len(a-C.O)-C.r)<=0;}//判断点A是否在圆C内
Circle MinCircle(Point *P,int n){//【求点集P的最小覆盖圆】
    mt19937 rd(time(0));shuffle(P+1,P+n+1,rd);Circle C;
    for(int i=2;i<=n;i++)if(!PIC(C,P[i]))
	{
        C=Circle(P[i],0);
        for(int j=1;j<i;j++)if(!PIC(C,P[j]))
		{
            C.O=(P[i]+P[j])*0.5,C.r=Len(P[j]-C.O);
            for(int k=1;k<j;k++)if(!PIC(C,P[k]))C=GetCircle(P[i],P[j],P[k]);
        }
    }
    return C;
}
const int N=1e5+3;
int n;
Point p[N];
int main()
{
	cin>>n;
	for(int i=1;i<=n;i++)p[i].In();
	Circle ans=MinCircle(p,n);
	cout<<fixed<<setprecision(10)<<ans.r<<endl;ans.O.Out();
}
posted @   Hanghang007  阅读(43)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· TypeScript + Deepseek 打造卜卦网站:技术与玄学的结合
· Manus的开源复刻OpenManus初探
· .NET Core 中如何实现缓存的预热?
· 三行代码完成国际化适配,妙~啊~
· 阿里巴巴 QwQ-32B真的超越了 DeepSeek R-1吗?
点击右上角即可分享
微信分享提示