计算几何模板

#include<algorithm>
#include<iostream>
#include<cstdlib>
#include<cstring>
#include<cstdio>
#include<string>
#include<cmath>
#include<ctime>
#include<queue>
#include<stack>
#include<map>
#include<set>
#define rre(i,r,l) for(int i=(r);i>=(l);i--)
#define re(i,l,r) for(int i=(l);i<=(r);i++)
#define Clear(a,b) memset(a,b,sizeof(a))
#define inout(x) printf("%d",(x))
#define douin(x) scanf("%lf",&x)
#define strin(x) scanf("%s",(x))
#define LLin(x) scanf("%lld",&x)
#define op operator
#define CSC main
typedef unsigned long long ULL;
typedef const int cint;
typedef long long LL;
using namespace std;
void inin(int &ret)
{
    ret=0;int f=0;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=1;ch=getchar();}
    while(ch>='0'&&ch<='9')ret*=10,ret+=ch-'0',ch=getchar();
    ret=f?-ret:ret;
}
const double eps=1e-8;
const double pi=acos(-1);
const double inf=2147483647;
double f(const long double &a){return a*a;}
inline int dcmp(const long double &a){return abs(a)<eps?0:a<0?-1:1;}
struct xl
{
    double x,y;
    xl(double x=0.,double y=0.):x(x),y(y){}
    void in(){douin(x),douin(y);}
    bool op == (const xl &rhs)const {return !dcmp(x-rhs.x)&&!dcmp(y-rhs.y);}
    bool op  < (const xl &rhs)const {return !dcmp(x-rhs.x)?y<rhs.y:x<rhs.x;}
    xl op + (const xl &rhs)const {return xl(x+rhs.x,y+rhs.y);}
    xl op - (const xl &rhs)const {return xl(x-rhs.x,y-rhs.y);}
    xl op * (const double &rhs)const {return xl(x*rhs,y*rhs);}
    xl op / (const double &rhs)const {return xl(x/rhs,y/rhs);}
    xl rotate(const double &rad)const {return xl(x*cos(rad)-y*sin(rad),x*sin(rad)+y*cos(rad));}
    double len()const {return sqrt(f(x)+f(y));}
    double angle()const {return atan2(y,x);}
};
double D_(const xl &a,const xl &b){return a.x*b.x+a.y*b.y;}
double X_(const xl &a,const xl &b){return a.x*b.y-a.y*b.x;}
double dis(const xl &a,const xl &b){return sqrt(f(a.x-b.x)+f(a.y-b.y));}
double dis2(const xl &a,const xl &b){return f(a.x-b.x)+f(a.y-b.y);}
double angle(const xl &a,const xl &b){return acos(D_(a,b)/a.len()/b.len());}
struct LINE
{
    xl u,v;double rad;
    LINE(xl u=xl(0.,0.),xl v=xl(0.,0.)):u(u),v(v){}
    void js(){rad=(v-u).angle();}
    bool op < (const LINE &rhs)const {return rad<rhs.rad;}
    double k()const {return !dcmp(u.x-v.x)?inf:(v.y-u.y)/(v.x-u.x);}
    double len()const {return (v-u).len();}
};
xl getjd(const LINE &a,const LINE &b)
{
    xl u=a.u-b.u,v1=a.v-a.u,v2=b.v-b.u;
    long double t=X_(v2,u)/X_(v1,v2);
    return a.u+v1*t;
}
bool pdxj(const LINE &a,const LINE &b)
{
    if(!dcmp(X_(a.u,b.v-b.u)-X_(a.v,b.v-b.u)))return 0;
    return 1;
}
double distoline(const xl &a,const LINE &l)
{
    return abs(X_(a-l.u,l.v-l.u))/l.len();
}
double distoseg(const xl &a,const LINE &l)
{
    xl v1=a-l.u,v2=a-l.v,u=l.v-l.u;
    int t1=dcmp(X_(v1,u)),t2=dcmp(X_(v2,u));
    if(t1>0&&t2>0)return v2.len();
    if(t1<0&&t2<0)return v1.len();
    return distoline(a,l);
}
xl getpoty(const xl &a,const LINE &l)
{
    xl v=l.v-l.u;
    return l.u+v*D_(a-l.u,v)/D_(v,v);
}
bool pdxjseg(const LINE &a,const LINE &b)
{
    xl u1=a.u-b.u,u2=a.v-b.u,u=a.v-a.u;
    xl v1=b.u-a.u,v2=b.v-a.u,v=b.v-b.u;
    int t1=dcmp(X_(u1,v)),t2=dcmp(X_(u2,v));
    int t3=dcmp(X_(v1,u)),t4=dcmp(X_(v2,v));
    return t1*t2<=0&&t3*t4<=0;
}
bool onseg(const xl &a,const LINE &l)
{
    return !dcmp(D_(a-l.u,a-l.v))&&dcmp(X_(a-l.u,a-l.v))<=0;
}
struct cir
{
    xl c;double r;
    cir(xl c=xl(0.,0.),double r=0):c(c),r(r){}
    bool op < (const cir &rhs)const {return c==rhs.c?r<rhs.r:c<rhs.c;}
    xl point (const double &rad)const {return c+xl(r*cos(rad),r*sin(rad));}
};
int getjd(const LINE &l,const cir &c,xl &ans1,xl &ans2)
{
    double d=distoline(c.c,l);
    if(dcmp(d-c.r)>0)return 0;
    if(!dcmp(d-c.r))return ans1=getpoty(c.c,l),1;
    xl temp=getpoty(c.c,l);
    xl v=c.c-temp;swap(v.x,v.y);
    double len=sqrt(f(c.r)-f(v.len()));
    v=v*len/v.len();
    ans1=temp+v,ans2=temp-v;
    return 2;
}
int getjd(const cir &a,const cir &b,xl &ans1,xl &ans2)
{
    double d=(b.c-a.c).len();
    if(!dcmp(d))if(!dcmp(b.r-a.r))return inf;else return 0;else ;
    if(dcmp(a.r+b.r-d)<0||dcmp(fabs(a.r-b.r)-d)>0)return 0;
    double rad=(b.c-a.c).angle();
    double da=acos((f(a.r)+f(d)-f(b.r))/(2*a.r*d));
    ans1=a.point(rad-da);
    xl ans3=a.point(rad+da);
    if(ans1==ans3)return 1;
    ans2=ans3;return 2;
}
int gettangents(cir &A,cir &B,LINE *ans)
{
    if(A.c==B.c)if(!dcmp(A.r-B.r))return inf;
    else return 0;else ;
    if(A.r<B.r)swap(A,B);
    double dd=f(A.c.x-B.c.x)+f(A.c.y-B.c.y);
    double rdiff=A.r-B.r,rsum=A.r+B.r;
    if(dcmp(dd-f(rdiff))<0)return 0;
    xl u=B.c-A.c;
    double base=atan2(B.c.y-A.c.y,B.c.x-A.c.x);
    if(!dcmp(dd-f(rdiff)))
    {
        xl p=A.point(base);
        ans[1]=LINE(p,p+u.rotate(pi/2));return 1;
    }
    double rad=acos((A.r-B.r)/sqrt(dd));
    ans[1]=LINE(A.point(base+rad),B.point(base+rad));
    ans[2]=LINE(A.point(base-rad),B.point(base-rad));
    if(dcmp(dd-f(rsum))<0)return 2;
    if(!dcmp(dd-f(rsum)))
    {
        xl p=A.point(base);
        ans[3]=LINE(p,p+u.rotate(pi/2));
        return 3;
    }
    rad=acos((A.r+B.r)/sqrt(dd));
    ans[3]=LINE(A.point(base+rad),B.point(base+rad+pi));
    ans[4]=LINE(A.point(base-rad),B.point(base-rad+pi));
    return 4;
}
cir outoftriangle(const xl &A,const xl &B,const xl &C)
{
    double bx=B.x-A.x,by=B.y-A.y;
    double cx=C.x-A.x,cy=C.y-A.y;
    double area=2*(bx*cy-by*cx);
    double x=(cy*(f(bx)+f(by))-by*(f(cx)+f(cy)))/area+A.x;
    double y=(bx*(f(cx)+f(cy))-cx*(f(bx)+f(by)))/area+A.y;
    return cir(xl(x,y),(A-xl(x,y)).len());
}
cir inoftriangle(const xl &A,const xl &B,const xl &C)
{
    double a=(B-C).len(),b=(C-A).len(),c=(A-B).len();
    xl cc=(A*a+B*b+C*c)/(a+b+c);
    return cir(cc,distoline(cc,LINE(A,B)));
}
struct pol
{
    xl a[100010];int n;
    void in(){inin(n);re(i,0,n-1)a[i].in();}
    void Sort(){sort(a,a+n);}
    xl& op [] (int x){return a[x];}
};
int pdinpol(const xl &p,pol &a)
{
    int sum=0,n=a.n;
    re(i,0,n)
    {
        if(onseg(p,LINE(a[i],a[(i+1)%n])))return -1;
        int k=dcmp(X_(a[(i+1)%n]-a[i],p-a[i]));
        int d1=dcmp(a[i].y-p.y);
        int d2=dcmp(a[(i+1)%n].y-p.y);
        if(k>0&&d1<=0&&d2>0)sum++;
        if(k<0&&d2<=0&&d1>0)sum--;
    }
    if(sum!=0)return 1;
    return 0;
}
void hull(pol &p,pol &ch)
{
    int n=p.n,m=0;
    p.Sort();
    re(i,0,n-1)
    {
        while(m>1&&dcmp(X_(ch[m-1]-ch[m-2],p[i]-ch[m-2]))<0)m--;
        ch[m++]=p[i];
    }
    int k=m;
    rre(i,n-1,0)
    {
        while(m>k&&dcmp(X_(ch[m-1]-ch[m-2],p[i]-ch[m-2]))<0)m--;
        ch[m++]=p[i];
    }
    if(n>1)m--;
    ch.n=m;
}
double lenofpol(pol &p)
{
    double ret=0;
    re(i,0,p.n-1)ret+=(p[(i+1)%p.n]-p[i]).len();
    return ret;
}
double areaofpol(pol &p)
{
    double ret=0;int n=p.n;
    re(i,1,n-2)ret+=X_(p[i]-p[0],p[i+1]-p[0]);
    return ret/2;
}
double diameterofhull(pol &p)
{
    int n=p.n;double ret=0.;
    if(n<=1)return 0.;
    if(n==2)return dis2(p[0],p[1]);
    int r=1;
    re(l,0,n-1)while(1)
    {
        double d=X_(p[l+1]-p[l],p[r+1]-p[r]);
        if(dcmp(d)<=0)
        {
            ret=max(ret,dis2(p[l],p[r]));
            if(!dcmp(d))ret=max(ret,dis2(p[l],p[r+1]));
            break;
        }
        r++,r%=n;
    }
    return ret;
}
bool onleft(const xl &a,const LINE &l){return X_(l.v-l.u,a-l.u)>0;}
bool halfplanej(LINE *l,int n,pol &ch)
{
    re(i,0,n-1)l[i].js();
    sort(l,l+n);
    int ll=0,rr=0;
    xl p[n];
    LINE q[n];q[0]=l[0];
    re(i,1,n-1)
    {
        while(ll<rr&&!onleft(p[rr-1],l[i]))rr--;
        while(ll<rr&&!onleft(p[ll],l[i]))ll++;
        q[++rr]=l[i];
        if(!dcmp(X_(q[rr].v-q[rr].u,q[rr-1].v-q[rr-1].u)))
        {
            rr--;
            if(onleft(l[i].u,q[rr]))q[rr]=l[i];
        }
        if(ll<rr)p[rr-1]=getjd(q[rr-1],q[rr]);
    }
    while(ll<rr&&!onleft(p[rr-1],q[ll]))rr--;
    if(rr-ll<=1)return 0;
    p[rr]=getjd(q[rr],q[ll]);
    re(i,ll,rr)ch[ch.n++]=p[i];
    ch[ch.n]=ch[0];return 1;
}
int main()
{
     return 0;
}