计算几何大模板

有错指出，长期更新

#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;
}