几何模板

#include <iostream>
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <vector>
#include <complex>
#include <queue>
#define eps 1e-10
#define pi acos(-1.0)
#define INF 0x7ffffff
#define maxn 50005
#define maxm 50005
#define next(i) ((i+1)%n)
using namespace std;
int cmp(double x)
{
    if(fabs(x)<eps)
        return 0;
    if(x>0)
        return 1;
    return -1;
}
inline double sqr(double x)
{
    return x*x;
}
struct point//点
{
    double x,y;
    point() {}
    point(double a,double b):x(a),y(b) {}
    void input()
    {
        scanf("%lf%lf",&x,&y);
    }
    friend point operator + (const point &a,const point &b)
    {
        return point(a.x+b.x,a.y+b.y);
    }
    friend point operator - (const point &a,const point &b)
    {
        return point(a.x-b.x,a.y-b.y);
    }
    friend bool operator == (const point &a,const point &b)
    {
        return cmp(a.x-b.x)==0 &&cmp(a.y-b.y)==0;
    }
    friend point operator * (const point &a,const double &b)
    {
        return point(a.x*b,a.y*b);
    }
    friend point operator * (const double &a,const point &b)
    {
        return point(a*b.x,a*b.y);
    }
    friend point operator / (const point &a,const double &b)
    {
        return point(a.x/b,a.y/b);
    }
    double norm()
    {
        return sqrt(sqr(x)+sqr(y));
    }
};
double det (const point &a,const point &b)
{
    return a.x*b.y-a.y*b.x;
}//叉积
double dot (const point &a,const point &b)
{
    return a.x*b.x+a.y*b.y;
}//点乘
double dist (const point &a,const point &b)
{
    return (a-b).norm();
}//距离
point rotate_point(const point &p,double A)
{
    double tx=p.x,ty=p.y;
    return point (tx*cos(A)-ty*sin(A),tx*sin(A)+ty*cos(A));
}//旋转
struct line
{
    point a,b;
    line() {}
    line(point x,point y):a(x),b(y) {}
};
line point_make_line(const point a,const point b)
{
    return line(a,b);
}//两个点a,b生成的一个线段
double dis_point_segment(const point p,const point s,const point t)
{
    if(cmp(dot(p-s,t-s))<0)
        return (p-s).norm();
    if(cmp(dot(p-t,s-t))<0)
        return (p-t).norm();
    return fabs(det(s-p,t-p)/dist(s,t));
}//p到线段st的距离
void PointProLine(const point p,const point s,const point t,point &cp)
{
    double r=dot((t-s),(p-s))/dot(t-s,t-s);
    cp=s+r*(t-s);
}//p到线段st的垂足，存在cp
bool PointOnSegment (point p,point s,point t)
{
    return cmp(det(p-s,t-s))==0&&cmp(dot(p-s,p-t))<=0;
}//点p是否在线段st上
bool parallel(line a,line b)
{
    return !cmp(det(a.a-a.b,b.a-b.b));
}//平行
bool line_make_point (line a,line b,point &res)
{
    if(parallel(a,b))
        return false;
    double s1=det(a.a-b.a,b.b-b.a);
    double s2=det(a.b-b.a,b.b-b.a);
    res=(s1*a.b-s2*a.a)/(s1-s2);
    return true;
}//交点
line move_d(line a,const double &len)
{
    point d=a.b-a.a;
    d=d/d.norm();
    d=rotate_point(d,pi/2);
    return line(a.a+d*len,a.b+d*len);
}
struct polygon//多边形
{
    int n;
    point a[maxn];
    polygon() {}
    double perimeter()
    {
        double sum=0;
        a[n]=a[0];
        for(int i=0; i<n; i++)
            sum+=(a[i+1]-a[i]).norm();
        return sum;
    }
    double area()
    {
        double sum=0;
        a[n]=a[0];
        for(int i=0; i<n; i++)
            sum+=det(a[i+1],a[i]);
        return sum/2.0;
    }
    int Point_In(point t)
    {
        int num=0,i,d1,d2,k;
        a[n]=a[0];
        for(int i=0; i<n; i++)
        {
            if(PointOnSegment(t,a[i],a[i+1]))
                return 2;
            k=cmp(det(a[i+1]-a[i],t-a[i]));
            d1=cmp(a[i].y=t.y);
            d2=cmp(a[i+1].y-t.y);
            if(k>0&&d1<=0&&d2>0)
                num++;
            if(k<0&&d2<=0&&d1>0)
                num--;
        }
        return num!=0;
    }//点是否在多边形内部
    point MassCenter()
    {
        point ans=point (0,0);
        if(cmp(area())==0)
            return ans;
        a[n]=a[0];
        for(int i=0; i<n; i++)
            ans=ans+(a[i]+a[i+1])*det(a[i+1],a[i]);
        return ans/area()/6.0;
    }//重心
    int Border_Int_Point_Num()
    {
        int num=0;
        a[n]=a[0];
        for(int i=0; i<n; i++)
            num+=__gcd (abs(int(a[i+1].x-a[i].x)),abs(int(a[i+1].y-a[i].y)));
        return num;
    }
    int Inside_Int_Point_Num()
    {
        return int(area())+1-Border_Int_Point_Num()/2;
    }//多边形内格点数
};
struct polygon_convex
{
    vector <point> P;
    polygon_convex(int size=0)
    {
        P.resize(size);
    }
};
bool comp_less(const point &a,const point &b)
{
    return cmp(a.x-b.x)<0||cmp(a.x-b.x)==0&&cmp(a.y-b.y)<0;
}
polygon_convex convex_hull(vector <point> a)
{
    polygon_convex res(2*a.size()+5);
    sort(a.begin(),a.end(),comp_less);
    a.erase(unique(a.begin(),a.end()),a.end());
    int m=0;
    for(int i=0; i<a.size(); i++)
    {
        while(m>1 && cmp(det(res.P[m-1]-res.P[m-2],a[i]-res.P[m-2]))<=0)
            m--;
        res.P[m++]=a[i];
    }
    int k=m;
    for(int i=int(a.size())-2; i>=0; i--)
    {
        while(m>k && cmp(det(res.P[m-1]-res.P[m-2],a[i]-res.P[m-2]))<=0)
            m--;
        res.P[m++]=a[i];
    }
    res.P.resize(m);
    if(a.size()>1)
        res.P.resize(m-1);
    return res;
}//凸包
bool contain(const polygon_convex &a,const point &b)
{
    int n=a.P.size();
    point g=(a.P[0]+a.P[n/3]+a.P[2*n/3])/3.0;
    int l=0,r=n;
    //二分凸包 g-a.P[a]-a.P[b]
    while(l+1<r)
    {
        int mid=(l+r)/2;
        if(cmp(det(a.P[l]-g,a.P[mid]-g))>0)
        {
            if(cmp(det(a.P[l]-g,b-g))>=0&&cmp(det(a.P[mid]-g,b-g))<0)
                r=mid;
            else l=mid;
        }
        else
        {
            if(cmp(det(a.P[l]-g,b-g))<0&&cmp(det(a.P[mid]-g,b-g))>=0)
                l=mid;
            else r=mid;
        }
    }
    r%=n;
    int z=cmp(det(a.P[r]-b,a.P[l]-b))-1;
    if(z==-2)
        return 1;
    return z;
}//O(longn)的判断点是否在凸包里
double convex_diameter (polygon_convex &a,int First,int Second)
{
    vector <point> p=a.P;
    int n=p.size();
    double maxd=0;
    if(n==1)
    {
        First=Second=0;
        return maxd;
    }
    for(int i=0,j=1; i<n; ++i)
    {
        while(cmp(det(p[next(i)]-p[i],p[j]-p[i])-det(p[next(i)]-p[i],p[next(j)]-p[i]))<0)
        {
            j=next(j);
        }
        double d=dist(p[i],p[j]);
        if(d>maxd)
        {
            maxd=d;
            First=i;
            Second=j;
        }
        d=dist(p[next(i)],p[next(j)]);
        if(d>maxd)
        {
            maxd=d;
            First=i;
            Second=j;
        }

    }
    return maxd;
}//凸包直径
inline bool operator < (const point &a,const point &b)
{
    return a.y+eps<b.y ||fabs(a.y-b.y)<eps &&a.x+eps<b.x;
}
inline double getA(const point &a)
{
    double res=atan2(a.y,a.x);
    if(res<0)
        res+=2*pi;
    return res;
}
point p[maxm],hull[maxm];
int n;
double w[maxm],sum[maxm];
inline void getHull()//逆时针凸包
{
    sort(p+1,p+n+1);
    int N=0;
    hull[++N]=p[1];
    for(int i=2; i<=n; i++)
    {
        while(N>1 &&det(hull[N]-hull[N-1],p[i]-hull[N-1])<=0)
            N--;
        hull[++N]=p[i];
    }
    int bak=N;
    for(int i=n-1; i>=1; i--)
    {
        while(N>bak&&det(hull[N]-hull[N-1],p[i]-hull[N-1])<=0)
            N--;
        hull[++N]=p[i];
    }
    n=N-1;
    for(int i=1; i<=n; i++)
        p[i+n]=p[i]=hull[i];
    p[n+n+1]=p[1];
    for(int i=1; i<=n; i++)
        w[i+n]=w[i]=getA(p[i+1]-p[i]);

    sum[0]=0;
    for(int i=1; i<=n*2; i++)
        sum[i]=sum[i-1]+det(p[i],p[i+1]);//预处理有向面积前缀和
}
inline int Find(double x)
{
    if(x<=w[1]||x>=w[n])
        return 1;
    return (upper_bound(w+1,w+n+1,x)-(w+1))+1;
}//找到第一个角度>x的边
point P,Q;
inline int getInter(int l,int r)
{
    int sign;
    if(det(Q-P,p[l]-P)<0)
        sign=-1;
    else sign=1;
    while(l+1<r)
    {
        int mid=(l+r)/2;
        if(det(Q-P,p[mid]-P)*sign>0)
            l=mid;
        else r=mid;
    }
    return r;
}//找到第一个和p[l]不在P->Q向量的同侧的点
inline point Intersect (const point &a,const point &b,const point &c,const point &d)
{
    double s1=det(c-a,b-a);
    double s2=det(d-a,b-a);
    return (c*s2-d*s1)/(s2-s1);
}
inline bool solve(point P,point Q)
{
    int i=Find (getA(Q-P));
    int j=Find (getA(P-Q));
    //两侧各找一点
    if(det(Q-P,p[i]-P)*det(Q-P,p[j]-P)>=0)//无交点
        return false;
    return true;
}
struct halfPlane
{
    double a,b,c;//ax+by+c<=0
    halfPlane(point p,point q)
    {
        a=p.y-q.y;
        b=q.x-p.x;
        c=det(p,q);
    }
    halfPlane(double aa,double bb,double cc)
    {
        a=aa;
        b=bb;
        c=cc;
    }
};
//计算点a和b连线与半平面L的交点
double calc(halfPlane &L,point &a)
{
    return a.x*L.a+a.y*L.b+L.c;
}
point Intersect(point &a,point &b,halfPlane &L)
{
    point res;
    double t1=calc(L,a),t2=calc(L,b);
    res.x=(t2*a.x-t1*b.x)/(t2-t1);
    res.y=(t2*a.y-t1*b.y)/(t2-t1);
    return res;
}
//将一个凸多边形和一个半平面交
polygon_convex cut(polygon_convex &a,halfPlane &L)
{
    int n=a.P.size();
    polygon_convex res;
    for(int i=0; i<n; i++)
    {
        if(calc(L,a.P[i])<-eps)
            res.P.push_back(a.P[i]);
        else
        {
            int j;
            j=i-1;
            if(j<0)
                j=n-1;
            if(calc(L,a.P[j])<-eps)
                res.P.push_back(Intersect(a.P[j],a.P[i],L));
            j=i+1;
            if(j==n)
                j=0;
            if(calc(L,a.P[j])<-eps)
                res.P.push_back(Intersect(a.P[i],a.P[j],L));
        }
    }
    return res;
}
polygon_convex core(polygon &a)
{
    polygon_convex res;
    res.P.push_back(point (-INF,-INF));
    res.P.push_back(point(INF,-INF));
    res.P.push_back(point(INF,INF));
    res.P.push_back(point(-INF,INF));
    int n=a.n;
    for(int i=0;i<n;i++)
    {
        halfPlane L(a.a[i],a.a[(i+1)%n]);
        res=cut(res,L);
    }
    return res;
}
typedef complex < double > Point;
typedef pair < Point , Point > Halfplane;
typedef vector< Point > Convex;
inline double dist(Halfplane hal)
{
    return sqrt((hal.second.imag()-hal.first.imag())*(hal.second.imag()-hal.first.imag())+(hal.second.real()-hal.first.real())*(hal.second.real()-hal.first.real()));
}
inline int sgn(double n)
{
    return fabs(n)<eps?0:(n<0?-1:1);
}
inline double cross(Point a,Point b)
{
    return (conj(a)*b).imag();
}
inline double dot(Point a,Point b)
{
    return(conj(a)*b).real();
}
inline double satisfy(Point a,Halfplane p)
{
    return sgn(cross(a-p.first,p.second-p.first))<=0;
}
Point crosspoint (const Halfplane &a,const Halfplane &b)
{
    double k=cross(b.first-b.second,a.first-b.second);
    k=k/(k-cross(b.first-b.second,a.second-b.second));
    return a.first+(a.second-a.first)*k;
}
bool cmp2(const Halfplane &a,const Halfplane &b)
{
    int res= sgn(arg(a.second-a.first)-arg(b.second-b.first));
    return res==0?satisfy(a.first,b):res<0;
}
vector < Point > halfplaneIntersection(vector<Halfplane>v)
{
    sort(v.begin(),v.end(),cmp2);
    deque < Halfplane > q;
    deque < Point > ans;
    q.push_back(v[0]);
    for(int i=1; i<int(v.size()); i++)
    {
        if(sgn(arg(v[i].second - v[i].first)-arg(v[i-1].second - v[i-1].first))==0)
            continue;
        while(ans.size()>0 && !satisfy(ans.back(),v[i]))
        {
            ans.pop_back();
            q.pop_back();
        }
        while(ans.size()>0 && !satisfy(ans.front(),v[i]))
        {
            ans.pop_front();
            q.pop_front();
        }
        ans.push_back(crosspoint (q.back(),v[i]));
        q.push_back(v[i]);
    }
    while(ans.size() > 0 && !satisfy(ans.back(),q.front()))
    {
        ans.pop_back();
        q.pop_back();
    }
    while(ans.size() > 0 && !satisfy(ans.front(),q.back()))
    {
        ans.pop_front();
        q.pop_front();
    }
    ans.push_back(crosspoint(q.back(),q.front()));
    return vector < Point > (ans.begin(),ans.end());
}
Convex convexIntersection(Convex v1,Convex v2)
{
    vector <Halfplane> h;
    for(int i=0; i<int(v1.size()); i++)
        h.push_back(Halfplane(v1[i],v1[(i+1) % v1.size()]));
    for(int i=0; i<int(v2.size()); i++)
        h.push_back(Halfplane(v2[i],v2[(i+1) % v2.size()]));
    return halfplaneIntersection(h);
}//凸多边形交
int dcmp(double k)
{
    return k<-eps?-1:k>eps?1:0;
}
double abs(const point &o)
{
    return sqrt(dot(o,o));
}
point crosspt(const point &a,const point &b,const point &p,const point &q)
{
    double a1=det(b-a,p-a);
    double a2=det(b-a,q-a);
    return (p*a2-q*a1)/(a2-a1);
}
double mysqrt(double n)
{
    return sqrt(max(0.0,n));
}
void circle_cross_line(point a,point b,point o,double r,point ret[],int &num)
{
    double x0=o.x,y0=o.y;
    double x1=a.x,y1=a.y;
    double x2=b.x,y2=b.y;
    double dx=x2-x1,dy=y2-y1;
    double A=dx*dx+dy*dy;
    double B=2*dx*(x1-x0) +2*dy*(y1-y0);
    double C=sqr(x1-x0)+sqr(y1-y0)-sqr(r);
    double delta =B*B-4*A*C;
    num=0;
    if(dcmp(delta) >= 0)
    {
        double t1=(-B-mysqrt(delta))/(2*A);
        double t2=(-B+mysqrt(delta))/(2*A);
        if(dcmp(t1-1)<=0&&dcmp(t1)>=0)//如果是直线，去掉这行if
        ret[num++]=point(x1+t1*dx,y1+t1*dy);
        if(dcmp(t2-1)<=0&&dcmp(t2)>=0)//如果是直线，去掉这行if
        ret[num++]=point(x1+t2*dx,y1+t2*dy);
    }
}//圆和直线交，ab(线段或直线），o圆心，r半径，ret交点，num交点数
point res[maxn];
double r;
double sector_area(const point &a , const point &b)
{
    double theta =atan2(a.y,a.x) -atan2(b.y,b.x);
    while(theta<=0)theta+=2*pi;
    while(theta>2*pi)theta-=2*pi;
    theta=min(theta,2*pi-theta);
    return r*r*theta/2;
}
double calc1(const point &a,const point &b)
{
    point p[2];
    int num=0;
    int ina=dcmp(abs(a)-r)<0;
    int inb=dcmp(abs(b)-r)<0;
    if(ina)
    {
        if(inb)
            return fabs(det(a,b))/2.0;
        else
        {
            circle_cross_line(a,b,point(0,0),r,p,num);
            return sector_area(b,p[0])+fabs(det(a,p[0]))/2.0;
        }
    }
    else
    {
        if(inb)
        {
            circle_cross_line(a,b,point(0,0),r,p,num);
            return sector_area(p[0],a)+fabs(det(p[0],b))/2.0;
        }
        else
        {
            circle_cross_line(a,b,point(0,0),r,p,num);
            if(num==2)
            {
                return sector_area(a,p[0])+sector_area(p[1],b)+fabs(det(p[0],p[1]))/2.0;
            }
            else return sector_area(a,b);
        }
    }
}
double area()
{
    double ret=0;
    for(int i=0;i<n;i++)
    {
        int sgn=dcmp(det(res[i],res[i+1]));
        if(sgn!=0)
        ret+=sgn*calc1(res[i],res[i+1]);
    }
    return ret;
}//res[n]=res[0],n全局变量，多边形点数，res，全局变量，逆时针存入多边形的所有点，r，全局变量，圆的半径