欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

计算几何大模板

程序员文章站 2024-01-14 17:35:28
...

有错指出,长期更新


#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<complex>
#include<vector>
#define db double
using namespace std;

namespace Complex{
    typedef complex<double> Point;
    typedef Point Vector;
    db dot  (Vector a,Vector b){return real(conj(a)*b);}
    db cross(Vector a,Vector b){return imag(conj(a)*b);}
    Vector rotate(Vector a,db rad){return a*exp(Point(0,rad));}
}

namespace geometry{
db eps=1e-8;
int dcmp(db x){if(fabs(x)<eps)return 0;return x<0?-1:1;}
const db pi=acos(-1);

struct Point{
    db x,y;
    Point(db a=0,db b=0):x(a),y(b){}
    void in(){scanf("%lf%lf",&x,&y);}
};

typedef Point Vector;

Vector operator +(Vector a,Vector b){return Vector(a.x+b.x,a.y+b.y);}
Vector operator -(Point  a,Point  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);    }
Vector operator /(Vector a,db     b){return Vector(a.x/b,a.y/b);    }

bool operator < (Point a,Point b){return a.x<b.x||(a.x==b.x&&a.y<b.y);      }
bool operator ==(Point a,Point b){return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;}
/*求叉积和点积*/
db cross  (Vector a,Vector b){return a.x*b.y-a.y*b.x;}
db dot    (Vector a,Vector b){return a.x*b.x+a.y*b.y;}
/*求向量长度*/
db length (Vector a)         {return sqrt(dot(a,a)); }
/*求两个向量的夹角*/
db angle  (Vector a,Vector b){return acos(dot(a,b)/length(a)/length(b));}
/*求三角形面积*/
db area   (Point a,Point b,Point c){return cross(b-a,c-a);              }
db area2  (Point &a,Point &b,Point &c){return cross(b-a,c-a);           }
/*求法线*/
Vector normal(Vector a)      {db l=length(a);return Vector(-a.y/l,a.x/l);}
/*旋转向量*/
Vector rotate(Point a,db rad){return Vector(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));}
/*求直线交点*/
Point getlinecut(Point a,Vector u,Point b,Vector v){
    Vector w=a-b;
    db     t=cross(v,w)/cross(u,v);
    return a+u*t;
}
/*点到直线距离*/
db distoline(Point p,Point a,Point b){
    Vector v1=b-a,v2=p-a;
    return fabs(cross(v1,v2)/length(v1));
}
/*求点到线段距离*/
db distosegm(Point p,Point a,Point b){
    if(a==b) return length(p-a);
    Vector v1=b-a,v2=p-a,v3=p-b;
    if(dcmp(dot(v1,v2))<0) return length(v2);
    if(dcmp(dot(v1,v3))>0) return length(v3);
    return fabs(cross(v1,v2)/length(v1));
}
/*求点到直线的垂足*/
Point getlinepro(Point p,Point a,Point b){
    Vector v=b-a;
    return a+v*(dot(v,p-a)/dot(v,v));
}
/*判断线段是否规范相交*/
bool segmpropcut(Point a1,Point a2,Point b1,Point b2){
    db c1=cross(a2-a1,b1-a1),c2=cross(a2-a1,b2-a1);
    db c3=cross(b2-b1,a1-b1),c4=cross(b2-b1,a2-b1);
    return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;
}
/*判断点是否在线段上*/
bool onsegm(Point p,Point a,Point b){
    return dcmp(cross(a-p,b-p))==0&&dcmp(dot(a-p,b-p))<0;
}
/*求多边形的面积*/
db polyarea(Point *p,int n){
    db ans=0;
    for(int i=2;i<=n-1;i++)ans+=area(p[1],p[i],p[i+1]);
    return ans/2;
}
/*求极角*/
db pangle(Vector v){return atan2(v.y,v.x);}
/*直线*/
struct line{
    Vector v;
    Point p;
    db ang;
    line(){}
    line(Point a,Vector b):p(a),v(b){ang=atan2(v.y,v.x);}
    Point getp(db t){return p+v*t;}
    bool operator <(line a)const{return ang<a.ang;}
};
Point getlinecut(line a,line b){return getlinecut(a.p,a.v,b.p,b.v);}
/*线段*/
struct segment{
    Point s,t;
    Vector v;
    db ang;
    segment(){}
    segment(Point a,Point b):s(a),t(b){v=b-a;ang=atan2(v.y,v.x);}
};
bool segmpropcut(segment a,segment b){return segmpropcut(a.s,a.t,b.s,b.t);}
bool onsegm(Point p,segment seg){return onsegm(p,seg.s,seg.t);}
/*圆*/
struct circle{
    Point o;
    db r;
    circle(Point a,db b):o(a),r(b){}
    Point getp(db rad){return Point(o.x+cos(rad)*r,o.y+sin(rad)*r);}
};
/*求直线与圆的交点,ans中存交点,返回交点个数*/
int getlinecirclecut(line l,circle cir,vector <Point> &ans){
    db a=l.v.x,b=l.p.x-cir.o.x,c=l.v.y,d=l.p.y-cir.o.y;
    db e=a*a+c*c,f=2*(a*b+c*d),g=b*b+d*d-cir.r*cir.r;
    db delta=f*f-4*e*g;
    if(dcmp(delta)<0) return 0;
    if(dcmp(delta)==0) {
        db t1=-f/(2*e);
        ans.push_back(l.getp(t1));
        return 1;
    }
    db t1=(-f-sqrt(delta))/(2*e);ans.push_back(l.getp(t1));
    db t2=(-f+sqrt(delta))/(2*e);ans.push_back(l.getp(t2));
    return 2;
}
int getlinecirclecut(line l,circle cir,db &t1,db &t2,vector <Point> &ans){
    db a=l.v.x,b=l.p.x-cir.o.x,c=l.v.y,d=l.p.y-cir.o.y;
    db e=a*a+c*c,f=2*(a*b+c*d),g=b*b+d*d-cir.r*cir.r;
    db delta=f*f-4*e*g;
    if(dcmp(delta)<0) return 0;
    if(dcmp(delta)==0) {
        t1=t2=-f/(2*e);
        ans.push_back(l.getp(t1));
        return 1;
    }
    t1=(-f-sqrt(delta))/(2*e);ans.push_back(l.getp(t1));
    t2=(-f+sqrt(delta))/(2*e);ans.push_back(l.getp(t2));
    return 2;
}
/*求两圆相交的交点,答案存在ans里,返回交点个数,-1表示圆重合*/
int getcircut(circle a,circle b,vector <Point> &ans){
    db d=length(a.o-b.o);
    if(dcmp(d)==0){if(dcmp(a.r-b.r)==0) return -1;return 0;}
    if(dcmp(a.r+b.r-d)<0) return 0;
    if(dcmp(fabs(a.r-b.r)-d)>0) return 0;
    db rad=pangle(a.o-b.o);
    db drad=acos((a.r*a.r+d*d-b.r*b.r)/(2*a.r*d));
    Point p1=a.getp(rad-drad),p2=a.getp(rad+drad);
    ans.push_back(p1);
    if(p1==p2) return 1;
    ans.push_back(p2);
    return 2;
}
/*求过一点到圆的切线,返回切线条数,存在ans里*/
int gettan(Point a,circle c,Vector *ans){
    Vector u=c.o-a;
    db dist=length(u);
    if(dist<c.r) return 0;
    if(dcmp(dist-c.r)==0){
        ans[0]=rotate(u,pi/2);
        return 1;
    }
    db ang=asin(c.r/dist);
    ans[0]=rotate(u,-ang);
    ans[1]=rotate(u,+ang);
    return 2;
}
int gettan(Point a,circle c,vector <line> &ans){
    Vector u=c.o-a;
    db dist=length(u);
    if(dist<c.r) return 0;
    if(dcmp(dist-c.r)==0){
        ans.push_back(line(a,rotate(u,pi/2)));
        return 1;
    }
    db ang=asin(c.r/dist);
    ans.push_back(line(a,rotate(u,-ang)));
    ans.push_back(line(a,rotate(u,+ang)));
    return 2;
}
db dis2(Point a,Point b){return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);}
/*求两个圆的公切线,-1为无穷条,a[i],b[i]分别为两圆上第i条切线的切点*/
int getcirtan(circle c1,circle c2,Point *a,Point *b){
    int cnt=0;
    if(c1.r<c2.r) swap(c1,c2),swap(a,b);
    Point o1=c1.o,o2=c2.o;
    db d2=dis2(o1,o2);
    db rdiff=c1.r-c2.r;
    if(d2<rdiff*rdiff) return 0;
    db base=pangle(o2-o1);
    if(d2==0&&c1.r==c2.r) return -1;
    if(d2==rdiff*rdiff) {
        a[++cnt]=c1.getp(base);b[cnt]=c2.getp(base);return 1;
    }
    db rsum=c1.r+c2.r;
    db ang=acos((c1.r-c2.r)/sqrt(d2));
    a[++cnt]=c1.getp(base+ang);b[cnt]=c2.getp(base+ang);
    a[++cnt]=c1.getp(base-ang);b[cnt]=c2.getp(base-ang);
    if(d2==rsum*rsum){
        ang=acos((c1.r+c2.r)/sqrt(d2));
        a[++cnt]=c1.getp(base+ang);b[cnt]=c2.getp(base+ang);
        a[++cnt]=c1.getp(base-ang);b[cnt]=c2.getp(base-ang);
    }
    return cnt;
}
int getcirtan(circle a,circle b,vector <line> &ans){
    Point *a1=new Point[5];
    Point *b1=new Point[5];
    int now=getcirtan(a,b,a1,b1);
    for(int i=1;i<=now;i++){ans.push_back(line(a1[i],a1[i]-b1[i]));}
    delete [] a1;
    delete [] b1;
    return now;
}
/*求三角形的外接圆*/
circle outcircle(Point p1,Point p2,Point p3){
  db Bx = p2.x-p1.x, By = p2.y-p1.y;
  db Cx = p3.x-p1.x, Cy = p3.y-p1.y;
  db D = 2*(Bx*Cy-By*Cx);
  db cx = (Cy*(Bx*Bx+By*By) - By*(Cx*Cx+Cy*Cy))/D + p1.x;
  db cy = (Bx*(Cx*Cx+Cy*Cy) - Cx*(Bx*Bx+By*By))/D + p1.y;
  Point p = Point(cx, cy);
  return circle(p, length(p1-p));
}
/*求三角形的内切圆*/
circle inscircle(Point p1,Point p2,Point p3){
  db a = length(p2-p3);
  db b = length(p3-p1);
  db c = length(p1-p2);
  Point p = (p1*a+p2*b+p3*c)/(a+b+c);
  return circle(p, distoline(p, p1, p2));
}
/*多边形*/
struct Poly{
    Point *p;
    db are;
    int sze,newsze;
    bool isare;
    void resize(){
        if(sze<newsze-1) return;
        Point *now=p;
        newsze=newsze<<1|1;
        p=new Point[newsze];
        for(int i=1;i<=sze;i++)p[i]=now[i];
        delete [] now;
    }
    Poly(){isare=0;}
    Poly(int size){p=new Point[size];newsze=size;isare=0;}
    void insert(Point a){
        resize();
        p[++sze]=a;isare=0;
    }
    int size(){return sze;}
    db getarea(){
        if(isare) return are;
        are=0;
        for(int i=2;i<=sze-1;i++){
            are+=area(p[1],p[i],p[i+1]);
        }
        are/=2;
        return are;
    }
    Point &operator [](int a){return p[a];}
    void clear(){sze=0;delete [] p;p=new Point[5];newsze=5;}
    ~Poly(){delete [] p;}
};
/*下标均从1开始*/
struct VPoly{
    vector <Point> p;
    void insert(Point a){p.push_back(a);}
    Point &operator [](int a){return p[a-1];}
    int size(){return p.size();}
    void clear(){p.clear();}
};
/*判断点与多边形关系,1内部,0外部,-1边上*/
int ispointinpoly(Point p,Poly poly){
    int wn=0;
    int n=poly.size();
    for(int i=1;i<=n;i++){
        if(onsegm(p,poly[i],poly[(i+1)%n])) return -1;
        int k =dcmp(cross(poly[(i+1)%n]-poly[i],p-poly[i]));
        int d1=dcmp(poly[i].y-p.y);
        int d2=dcmp(poly[(i+1)%n].y-p.y);
        if(k>0&&d1<=0&&d2>0) wn++;
        if(k>0&&d2<=0&&d1>0) wn--;
    }
    if(wn!=0) return 1;
    return 0;
}
/*求凸包,答案装在ans,若不希望凸包边上有点就用下一个函数*/
int convexhull(Point *p,int n,Point *ans){
    sort(p+1,p+n+1);
    int m=0;
    for(int i=1;i<=n;i++){
        while(m>2&&dcmp(cross(ans[m-1]-ans[m-2],p[i]-ans[m-2]))<=0) m--;
        ans[++m]=p[i];
    }
    int k=m;
    for(int i=n-1;i>=1;i--){
        while(m>=k&&dcmp(cross(ans[m-1]-ans[m-2],p[i]-ans[m-2]))<=0) m--;
        ans[++m]=p[i];
    }
    return m;
}
int convexhull_noside(Point *p,int n,Point *ans){
    sort(p+1,p+n+1);
    int m=0;
    for(int i=1;i<=n;i++){
        while(m>2&&dcmp(cross(ans[m-1]-ans[m-2],p[i]-ans[m-2]))<0) m--;
        ans[++m]=p[i];
    }
    int k=m;
    for(int i=n-1;i>=1;i--){
        while(m>=k&&dcmp(cross(ans[m-1]-ans[m-2],p[i]-ans[m-2]))<0) m--;
        ans[++m]=p[i];
    }
    return m;
}
/*下标从0开始*/
int convexhull_zero(Point *p,int n,Point *ans){
    sort(p,p+n);
    int m=0;
    for(int i=0;i<n;i++){
        while(m>1&&dcmp(cross(ans[m-1]-ans[m-2],p[i]-ans[m-2]))<=0) m--;
        ans[m++]=p[i];
    }
    int k=m;
    for(int i=n-2;i>=0;i--){
        while(m>=k&&dcmp(cross(ans[m-1]-ans[m-2],p[i]-ans[m-2]))<=0) m--;
        ans[m++]=p[i];
    }
    if(n>1)m--;
    return m;
}
int convexhull_zero_noside(Point *p,int n,Point *ans){
    sort(p,p+n);
    int m=0;
    for(int i=0;i<n;i++){
        while(m>1&&dcmp(cross(ans[m-1]-ans[m-2],p[i]-ans[m-2]))<=0) m--;
        ans[m++]=p[i];
    }
    int k=m;
    for(int i=n-2;i>=0;i--){
        while(m>k&&dcmp(cross(ans[m-1]-ans[m-2],p[i]-ans[m-2]))<=0) m--;
        ans[m++]=p[i];
    }
    if(n>1)m--;
    return m;
}
/*半平面交*/
/*判断点在直线的坐边,线上不算*/
bool onleft(line l,Point p){
    return dcmp(cross(l.v,p-l.p))>0;
}
/*半平面交可以为点时使用*/
bool onleft2(line l,Point p){
    return dcmp(cross(l.v,p-l.p))>=0;
}
/*返回点数*/
/*下标为0*/
int halfcut(line *l,int n,Point *ans){
    sort(l,l+n);
    int fi,la;
    Point *p= new Point[n];
    line *q= new line[n];
    q[fi=la=0]=l[0];
    for(int i=1;i<n;i++){
        while(fi<la&&!onleft(l[i],p[la-1])) la--;
        while(fi<la&&!onleft(l[i],p[fi]))   fi++;
        q[++la]=l[i];
        if(fabs(cross(q[la].v,q[la-1].v))<eps){
            la--;
            if(onleft(q[la],l[i].p)) q[la]=l[i];
        }
        if(fi<la) p[la-1]=getlinecut(q[la-1],q[la]);
    }
    while(fi<la&&!onleft(q[fi],p[la-1]))la--;
    if(la-fi<=1) return 0;
    p[la]=getlinecut(q[la],q[fi]);
    int m=0;
    for(int i=fi;i<=la;i++) ans[m++]=p[i];
    delete [] p;
    delete [] q;
    return m;
}
/*下标为1,可能有错*/
int halfcut_one(line *l,int n,Point *ans){
    sort(l+1,l+n+1);
    int fi,la;
    Point *p= new Point[n];
    line *q= new line[n];
    q[fi=la=1]=l[1];
    for(int i=2;i<=n;i++){
        while(fi<la&&!onleft(l[i],p[la-1])) la--;
        while(fi<la&&!onleft(l[i],p[fi]))   fi++;
        q[++la]=l[i];
        if(fabs(cross(q[la].v,q[la].v))<eps){
            la--;
            if(onleft(q[la],l[i].p)) q[la]=l[i];
        }
        if(fi<la) p[la-1]=getlinecut(q[la-1],q[la]);
    }
    while(fi<la&&!onleft(q[fi],p[la-1]))la--;
    if(la-fi<=1) return 0;
    p[la]=getlinecut(q[la],q[fi]);
    int m=0;
    for(int i=fi;i<=la;i++) ans[++m]=p[i];
    delete [] p;
    delete [] q;
    return m;
}
/*旋转卡壳*/
db dsts(Point a,Point b,Point c,Point d){
    return min(min(distosegm(c,a,b),distosegm(d,a,b)),min(distosegm(a,c,d),distosegm(b,c,d)));
}
/*求两个凸多边形之间的最短距离, U17652 岛屿连接:https://www.luogu.org/record/show?rid=5244376*/
db get_two_poly_min_dis(Point *p1,Point *p2,int n,int m){
    int pos1=0,pos2=0;
    db ans=1e90;
    for(int i=0;i<n;i++){
        if(p1[i].y>=p1[pos1].y){
            if(p1[i].y>p1[pos1].y) pos1=i;
            else if(p1[i].x>p1[pos1].x) pos1=i;
        }
    }
    for(int i=0;i<m;i++){
        if(p2[i].y<=p2[pos2].y){
            if(p2[i].y<p2[pos2].y) pos2=i;
            else if(p2[i].x<p2[pos2].x) pos2=i;
        }
    }
    db ls;
    for(int i=0;i<n;i++){
        while(ls=(cross(p1[(pos1+1)%n]-p1[pos1],p2[(pos2+1)%m]-p1[pos1])-cross(p1[(pos1+1)%n]-p1[pos1],p2[pos2]-p1[pos1]))>eps)
        pos2=(pos2+1)%m;
        if(ls<-eps) ans=min(ans,distosegm(p2[pos2],p1[(pos1+1)%n],p1[pos1]));
        else ans=min(ans,dsts(p1[(pos1+1)%n],p1[pos1],p2[(pos2+1)%m],p2[pos2]));
        pos1=(pos1+1)%n;
    }
    return ans;
}
db dis(Point a,Point b){return length(a-b);}
/*求凸多边形的直径*/
db RotatingCaliper(int m,Point *ch)
{
    int q = 1;
    ch[m] = ch[0];
    db ans = 0;
    for(int p = 0; p < m; p++)
    {
        while(fabs(cross(ch[q+1]-ch[p+1], ch[p]-ch[p+1])) > fabs(cross(ch[q]-ch[p+1], ch[p]-ch[p+1]))) q = (q+1)%m;
        ans = max(ans,max(dis(ch[p],ch[q]),dis(ch[p+1],ch[q+1])));
        ans = max(ans,max(dis(ch[p],ch[q+1]),dis(ch[p+1],ch[q])));
    }
    return ans;
}
/*角度转弧度*/
db torad(db deg){return deg/180*pi;}
/*******三维几何********/
/*经纬度转换为空间坐标*/
void get_coord(db r,db lat,db lng,db &x,db &y,db &z){
    lat=torad(lat);
    lng=torad(lng);
    x=r*cos(lat)*cos(lng);
    y=r*cos(lat)*sin(lng);
    z=r*sin(lat);
}
struct Point3{
    db x,y,z;
    Point3(db a=0,db b=0,db c=0):x(a),y(b),z(c){}
};
typedef Point3 Vector3;
Vector3 operator +(Vector3 a,Vector3 b){return Vector3(a.x+b.x,a.y+b.y,a.z+b.z);}
Vector3 operator -(Point3  a,Point3  b){return Vector3(a.x-b.x,a.y-b.y,a.z-b.z);}
Vector3 operator *(Vector3 a,db      b){return Vector3(a.x*b,  a.y*b,  a.z*b)  ;}
Vector3 operator /(Vector3 a,db      b){return Vector3(a.x/b,  a.y/b,  a.z/b)  ;}
db dot3(Vector3 a,Vector3 b){return a.x*b.x+a.y*b.y+a.z*b.z;}
Vector3 cross(Vector3 a,Vector3 b){return Vector3(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x);}

}
using namespace geometry;
int main(){
    return 0;
}