计算几何

没有封完的全家桶

#include<bits/stdc++.h>
using namespace std;
#define ld double
#define vec Point
const ld eps=1e-10;
const int N=10;
int dcmp(ld a){return fabs(a)<eps?0:(a<0?-1:1);}
struct Point{
    ld x,y;Point(ld a=0,ld b=0){x=a,y=b;}
    friend ld operator*(const vec &a,const vec &b){return a.x*b.y-a.y*b.x;}
    friend ld operator^(const vec &a,const vec &b){return a.x*b.x+a.y*b.y;}
    friend vec operator+(vec a,vec b){return vec(a.x+b.x,a.y+b.y);}
    friend vec operator-(vec a,vec b){return vec(a.x-b.x,a.y-b.y);}
    vec operator*(ld a){return vec(a*x,a*y);}
    vec operator/(ld a){return vec(x/a,y/a);}
    bool operator==(Point a){return !dcmp(x-a.x)&&!dcmp(y-a.y);}
};
namespace calc{
    ld Len(vec a){return sqrt(a^a);}
    Point cross_LL(Point a,Point b,Point c,Point d){return a+(b-a)*(((d-c)*(a-c))/((b-a)*(d-c)));}
    bool judge_SegLine_P(Point p,vec a,vec b){return dcmp((p-a)^(p-b))<=0&&!dcmp((p-a)*(p-b));}
    bool judge_StraintLine_P(Point p,vec a,vec b){return !dcmp((p-a)*(p-b));}
    bool judge_StraintSegLine(Point a,Point b,Point c,Point d){return judge_SegLine_P(cross_LL(a,b,c,d),c,d);}
    bool judge_SegSegLine(Point a,Point b,Point c,Point d){
        ld c1((b-a)*(c-a)),c2((b-a)*(d-a)),c3((d-c)*(a-c)),c4((d-c)*(b-c));
        return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;
    }
    ld dis_PL(Point p,vec a,vec b){
        if((dcmp((p-a)^(b-a))<0)||(a==b)) return Len(p-a); 
        if(dcmp((p-b)^(b-a))>0) return Len(p-b);
        return fabs((p-b)^(b-a))/Len(b-a);
    }
    Point FootPoint(Point p,vec a,vec b){
        ld len1=((p-a)^(b-a))/Len(b-a),len2=-1.0*((p-b)^(b-a))/Len(b-a);
        return a+(b-a)*(len1/(len1+len2));
    }
    Point Symmetry(Point p,vec a,vec b){return p+((FootPoint(p,a,b)-p)*2);}
    int pip(Point *p,int n,Point a){
        int cnt(0);ld tmp;
        for(int i=1;i<=n;++i){
            int j(i<n?i+1:1);
            if(judge_SegLine_P(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;
        }
        return cnt&1;
    }
    int judge_PoPo(Point *A,int n,Point *B,int m){
        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(judge_SegSegLine(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;
    }
	ld Po_Area(Point *p,int n){
		ld S(0);
		for(int i=1;i<=n;++i) S+=p[i]*p[i<n?i+1:1];
		return S/2.0;
	}
	int ConvexHull(Point *P,int n,Point *cp){
    	sort(P+1,P+n+1,[](vec a,vec b){return a.x==b.x?a.y<b.y:a.x<b.x;});int t=0;
    	for(int i=1;i<=n;++i) while(t>1&&dcmp((cp[t]-cp[t-1])*(P[i]-cp[t-1]))<=0) --t,cp[++t]=P[i];
		int St=t;
		for(int i=n-1;i>=1;--i) while(t>St&&dcmp((cp[t]-cp[t-1])*(P[i]-cp[t-1]))<=0) --t,cp[++t]=P[i];
		return --t;
	}
	vec V1[N],V2[N];
	int Mincowski(Point *P1,int n,Point *P2,int m,vec *V){//【闵可夫斯基和】求两个凸包{P1},{P2}的向量集合{V}={P1+P2}构成的凸包
		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]+(dcmp((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;
	}
}
namespace modify{
    Point turn_pp(Point a,Point b,ld theta){
        ld x((a.x-b.x)*cos(theta)+(a.y-b.y)*sin(theta)+b.x);
        ld y(-(a.x-b.x)*sin(theta)+(a.y-b.y)*cos(theta)+b.y);
        return Point(x,y);
    }
}