模板 - 计算几何

这个模板的缺点：有很多模板都是非常C语言化的，虽然可读性比较差但是应该性能非常感人。鉴于ACM可以开O2所以方向用STL就好了。但是有的地方不好，比如半平面交的排序，对同一个向量多次判断倾斜角，其实预处理的时候要是需要，就把倾斜角初始化就好了。

基础

//不要输出-0.0之类的数

const double eps=1e-8;
const double pi=acos(-1.0);

//判断浮点数的符号
inline int cmp(double x){
    return (fabs(x)<eps)?0:((x>0.0)?1:-1);
}

inline double sqr(double x){
    return x*x;
}

多边形

点

struct Point {
    double x,y;
    Point() {};
    Point(const double x,const double y):x(x),y(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 Point operator*(const Point &p,const double k) {
        return Point(p.x*k,p.y*k);
    }
    friend Point operator*(const double k,const Point &p) {
        return Point(p.x*k,p.y*k);
    }
    friend Point operator/(const Point &p,const double k) {
        return Point(p.x/k,p.y/k);
    }
    friend bool operator==(const Point &a,const Point &b) {
        return cmp(a.x-b.x)==0&&cmp(a.y-b.y)==0;
    }
    Point rotate(double A) {
        //向量绕原点旋转A弧度
        return Point(x*cos(A)-y*sin(A),x*sin(A)+y*cos(A));
    }
    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 sqrt(sqr(a.x-b.x)+sqr(a.y-b.y));
}

线段

struct Line {
    Point a,b;
    Line() {};
    Line(const Point &a,const Point &b):a(a),b(b) {};
    Line move_dist(const double &d) {
        //向法向平移d单位长度
        //单位法向量n，从a指向b
        Point n=b-a;
        n=n/n.norm();
        //左旋90度
        n=n.rotate(pi/2.0);
        return Line(a+n*d,b+n*d);
    }
};

double dist_point_to_line(const Point &p,const Line &l) {
    Point a=l.a,b=l.b;
    //当a与b可以重合时，这里要加上下面的语句
    /*if(a==b)
        return a.dist(p);*/
    if(cmp(dot(p-a,b-a))<0)
        return dist(p,a);
    if(cmp(dot(p-b,a-b))<0)
        return dist(p,b);
    return fabs(det(a-p,b-p)/dist(a,b));
}
Point point_project_on_line(const Point &p,const Line &l)  {
    Point a=l.a,b=l.b;
    double r=dot(b-a,p-a)/dot(b-a,b-a);
    return a+(b-a)*r;
}
bool point_on_line(const Point &p,const Line &l)  {
    Point a=l.a,b=l.b;
    //这里的line是线段
    //第一个cmp意思是叉积等于0，意味着直线穿过该点
    //第二个cmp的<=意思是点在线段内（含端点），当改为<为点在线段内（不含端点）
    return cmp(det(p-a,b-a))==0&&cmp(dot(p-a,p-b))<=0;
}
bool parallel(const Line &tl,const Line &l)  {
    Point a=tl.a,b=tl.b;
    //叉积等于0，意味着向量平行
    return !cmp(det(a-b,l.a-l.b));
}
bool intersect(const Line &tl,const Line &l,Point &p) {
    Point a=tl.a,b=tl.b;
    //判断直线是否相交，相交则求出交点（不需要交点可以直接return）
    if(parallel(tl,l))
        return false;
    double s1=det(a-l.a,l.b-l.a);
    double s2=det(b-l.a,l.b-l.a);
    p=(b*s1-a*s2)/(s1-s2);
    return true;
}

多边形

const int MAXN=105;
struct Polygon {
    int n;
    Point a[MAXN];
    Polygon() {};
    double perimeter() {
        double sum=0.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.0;
        a[n]=a[0];
        for(int i=0; i<n; i++)
            sum+=det(a[i+1],a[i]);
        return sum/2.0;
    }
    Point masscenter(){
        Point ans(0.0,0.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_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_point_num(){
        return (int)area()+1-border_point_num()/2;
    }
};

int point_in_polygon(Point &p,Polygon &po) {
    Point *a=po.a;
    int n=po.n;
    int num=0,d1,d2,k;

    a[n]=a[0];
    for(int i=0; i<n; i++) {
        if(point_on_line(p,Line(a[i],a[i+1])))
            return 2;
        k=cmp(det(a[i+1]-a[i],p-a[i]));
        d1=cmp(a[i].y-p.y);
        d2=cmp(a[i+1].y-p.y);
        if(k>0&&d1<=0&&d2>0)
            num++;
        if(k<0&&d2<=0&&d1>0)
            num--;
    }
    return num!=0;
}

凸包

struct Polygon_Convex {
    vector<Point> P;
    Polygon_Convex(int Size=0) {
        P.resize(Size);
    }
    Polygon to_polygon() {
        //注意多边形的最大点数要够
        Polygon p;
        p.n=P.size();
        for(int i=0; i<p.n; i++) {
            p.a[i]=P[i];
        }
        return p;
    }
    double diameter(int &First,int &Second){
        //旋转卡壳求直径，O(n)
        vector<Point> &p=P;
        int n=P.size();
        double maxd=0.0;
        if(n==1){
            First=Second=0;
            return maxd;
        }
        #define next(i) ((i+1)%n)
        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;
            }
        }
        #undef next(i)
        return maxd;
    }
};

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-(a.size()>1));
    return res;
}

int point_in_polygon_convex(const Point &p,const Polygon_Convex &pc) {
    //0在外部，1在内部，2在边界上
    //包括边界
    int n=pc.P.size();
    const vector<Point> &P=pc.P;
    //找一个内部点
    Point g=(P[0]+P[n/3]+P[2*n/3])/3.0;
    int l=0,r=n;
    while(l+1<r){
        int mid=(l+r)>>1;
        if(cmp(det(P[l]-g,P[mid]-g))>0){
            if(cmp(det(P[l]-g,p-g))>=0&&cmp(det(P[mid]-g,p-g))<0)
                r=mid;
            else
                l=mid;
        }
        else{
            if(cmp(det(P[l]-g,p-g))<0&&cmp(det(P[mid]-g,p-g))>=0)
                l=mid;
            else
                r=mid;
        }
    }
    r%=n;
    int z=cmp(det(P[r]-p,P[l]-p));
    //z==0在边界上，三点共线
    //z==1在凸包外
    //z==-1在凸包内
    return (z+2)%3;
}

半平面

强烈建议加上四周围的inf边，这样传入的v既不会是空的，也不会有相反向量叉积做除法（转180度之前会遇到inf边）。
小心三线共点返回的半平面，有点东西。（好像可以被反向平行检测掉）

    vector<Halfplane> Hp;
    Hp.push_back(Halfplane(Point(-inf,-inf),Point(inf,-inf)));
    Hp.push_back(Halfplane(Point(inf,-inf),Point(inf,inf)));
    Hp.push_back(Halfplane(Point(inf,inf),Point(-inf,inf)));
    Hp.push_back(Halfplane(Point(-inf,inf),Point(-inf,-inf)));

struct Halfplane {
    //向量first->second的左侧
    Point first,second;
    Halfplane() {};
    Halfplane(Point p1,Point p2):first(p1),second(p2) {};
};

inline int satisfy(Point a,Halfplane p) {
    return cmp(det(a-p.first,p.second-p.first))<=0;
}

Point intersect_point(const Halfplane &a,const Halfplane &b) {
    double k=det(b.first-b.second,a.first-b.second);
    double t=det(b.first-b.second,a.second-b.second);
    //反向向量已经被半平面交制裁了，这个函数调用之前就要先保证不平行，否则后果自负
    //把边界也放进来了，所以反向向量至少转180度会遇到边界
    k=k/(k-t);
    return a.first+(a.second-a.first)*k;
}

inline bool compare(const Halfplane &a,const Halfplane &b) {
    int res=cmp((a.second-a.first).angle()-(b.second-b.first).angle());
    return res==0?satisfy(a.first,b):res<0;
}

inline bool parallel(const Halfplane &a,const Halfplane &b){
    Point pa=a.second-a.first;
    Point pb=b.second-b.first;
    return !cmp(det(pa,pb));
}

//半平面交，O(nlogn)
vector<Point> halfplane_intersection(vector<Halfplane> v) {
    //半平面把边界放进来，不可能是空的
    if(v.empty())
        return vector<Point>();
    sort(v.begin(),v.end(),compare);
    deque<Halfplane> q;
    deque<Point> ans;
    q.push_back(v[0]);

    int vs=v.size();
    for(int i=1; i<vs; ++i) {
        if(cmp((v[i].second-v[i].first).angle()-(v[i-1].second-v[i-1].first).angle())==0)
            continue;
        while(!ans.empty()&&!satisfy(ans.back(),v[i])) {
            ans.pop_back();
            q.pop_back();
        }
        while(!ans.empty()&&!satisfy(ans.front(),v[i])) {
            ans.pop_front();
            q.pop_front();
        }
        //前面已经去掉平行了，再平行就是反平行
        if(parallel(q.back(),v[i]))
            return vector<Point>();
        ans.push_back(intersect_point(q.back(),v[i]));
        q.push_back(v[i]);
    }
    while(!ans.empty()&&!satisfy(ans.back(),q.front())) {
        ans.pop_back();
        q.pop_back();
    }
    while(!ans.empty()&&!satisfy(ans.front(),q.back())) {
        ans.pop_front();
        q.pop_front();
    }
    ans.push_back(intersect_point(q.back(),q.front()));
    return vector<Point>(ans.begin(),ans.end());
}