计算几何板子

        这里放计算几何的板子。关于原理，推荐一篇很好的文章。
        长期更新。

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <map>
#include <queue>
#include <set>
#include <stack>
#include <vector>
using namespace std;
const double PI = acos(-1.0);  //π
const double EPS = 1e-8;       //精度控制
const double INF = 1e100;      //无穷大
int sgn(double s)              //判断浮点数的符号
{
    if (fabs(s) < EPS) return 0;
    if (s > 0) return 1;
    return -1;
}
int cmp(double x, double y)  //判断x和y的大小关系,-1表示x<y
{
    if (fabs(x - y) < EPS) return 0;
    if (x < y) return -1;
    return 1;
}
double r2d(double rad)  //弧度转角度
{
    return rad / PI * 180.0;
}
double d2r(double degree)  //角度转弧度
{
    return degree / 180.0 * PI;
}
struct Point  //点与向量
{
public:
    double x, y;
    Point(double a = 0, double b = 0)
        : x(a), y(b) {}
    void input()  //输入
    {
        scanf("%lf%lf", &x, &y);
    }
    void debug()  //输出
    {
        printf("(%lf,%lf)\n", x, y);
    }
    Point operator+(Point b)  //向量相加
    {
        return Point(x + b.x, y + b.y);
    }
    Point operator-(Point b)  //向量相减
    {
        return Point(x - b.x, y - b.y);
    }
    Point operator*(double b)  //向量数乘
    {
        return Point(x * b, y * b);
    }
    double operator*(Point b)  //向量内积
    {
        return x * b.x + y * b.y;
    }
    Point operator/(double b)  //除法
    {
        return Point(x / b, y / b);
    }
    bool operator<(const Point& b)  //大小比较
    {
        if (cmp(x, b.x) == 0) return y < b.y;
        return x < b.x;
    }
    bool operator==(const Point& b)  //判等
    {
        return cmp(x, b.x) == 0 && cmp(y, b.y) == 0;
    }
    double cross(Point b) const  //向量叉乘
    {
        return x * b.y - y * b.x;
    }
    double angle()  //极角,弧度(-PI, PI]
    {
        return atan2(y, x);
    }
    double length()  //模
    {
        return sqrt((*this) * (*this));
    }
    double length2()
    {
        return (*this) * (*this);
    }
    double angleTo(Point b)  //向量角度,弧度
    {
        return acos(((*this) * b) / length() / b.length());
    }
    double angleTo(Point a, Point b)  //从该点看其它两点的角度,弧度
    {
        return Point(a - *this).angleTo(b - *this);
    }
    double area(Point a, Point b)  //求三点组成平行四边形有向面积
    {
        return (a - *this).cross(b - *this);
    }
    double dist(Point b)  //求两点间欧式距离
    {
        return (*this - b).length();
    }
    double dist2(Point b)
    {
        return (*this - b).length2();
    }
    double manhattanDis(Point b)  //求两点间曼哈顿距离
    {
        return fabs(x - b.x) + fabs(y - b.y);
    }
    Point rotate(Point p, double ang)  //绕点p逆时针旋转ang弧度,不修改本身的值
    {
        Point v = (*this) - p;
        double c = cos(ang), s = sin(ang);
        return p + Point(v.x * c - v.y * s, v.x * s + v.y * c);
    }
    Point rotate(double ang)  //直接逆时针旋转ang弧度
    {
        return rotate(Point(0, 0), ang);
    }
    Point rotLeft()  //左转90
    {
        return Point(-y, x);
    }
    Point rotRight()  //右转90
    {
        return Point(y, -x);
    }
    Point regular()  //化为单位向量
    {
        return *this / length();
    }
    Point normal()  //左转90的单位法向量
    {
        return rotLeft().regular();
    }
    Point trunc(double r)  //化为长度为r的向量
    {
        return regular() * r;
    }
    bool parrelTo(Point b)  //向量平行
    {
        return sgn(cross(b)) == 0;
    }
    bool verticalTo(Point b)  //向量垂直
    {
        return sgn((*this) * b) == 0;
    }
};
typedef Point Vector;
int toLeftTest(Point from, Point to, Point test)  //测试test在不在from->to的左边,1在左边,-1在右边,0共线
{
    return sgn((to - from).cross(test - to));
}
struct Circle  //圆
{
    Point o;   //圆心
    double r;  //半径
    Circle(Point o = Point(), double r = 0.0)
        : o(o), r(r) {}

    void input()
    {
        scanf("%lf%lf%lf", &o.x, &o.y, &r);
    }
    void debug()
    {
        printf("(%lf,%lf,%lf)\n", o.x, o.y, r);
    }
    Circle invertToCicle(Point P, double R)  //外部反演点反演,反演为圆; p为反演中心,r为反演半径
    {
        Circle ans;
        double d1 = P.dist(o);
        ans.r = r / (d1 * d1 - r * r) * R * R;
        double d2 = R * R / (d1 - r) - ans.r;
        ans.o = P + (o - P) * (d2 / d1);
        return ans;
    }
};
struct Line
{
    Point pos;
    Vector to;
    Line(Point p = Point(), Vector t = Point())  //顶点和方向向量
        : pos(p), to(t)
    {
    }
    static Line tp(Point a, Point b)  //两点构造
    {
        return Line(a, b - a);
    }
    Line(double x1, double y1, double x2, double y2)  //两点构造
    {
        pos = Point(x1, y1), to = Point(x2, y2) - Point(x1, y1);
    }
    Point point(double s)  //获取某处的一个点
    {
        return pos + to * s;
    }
    double angle()  //获得直线倾角,[0,PI)
    {
        return acos(to.regular().x * sgn(to.regular().y));
    }
    int onLine(Point s)  //判断点与直线的位置关系0在直线外;1在直线上
    {
        return sgn(to.cross(s - pos)) == 0;
    }
    bool operator==(Line l)  //两条直线是否是同一条(重合)
    {
        return l.onLine(point(0)) && l.onLine(point(1));
    }
    bool operator!=(Line l)
    {
        return !((*this) == l);
    }
    bool parrelTo(Line l)  //直线平行
    {
        return to.parrelTo(l.to);
    }
    bool verticalTo(Line l)  //直线垂直
    {
        return to.verticalTo(l.to);
    }
    Point insWithLine(Line l)  //与另一条直线的交点
    {
        return point(l.to.cross(pos - l.pos) / to.cross(l.to));
    }
    double disToPoint(Point s)  //点到直线的距离
    {
        return fabs(to.cross(s - pos) / to.length());
    }
    Point projection(Point s)  //点在直线上的投影
    {
        return point((to * (s - pos)) / (to * to));
    }
    Circle invertToCircle(Point P, double R)  //直线外一点反演为圆
    {
        Line lp(P, to.rotLeft());
        Point ps = this->insWithLine(lp);
        Circle c;
        c.r = R * R / this->disToPoint(P) / 2.0;
        c.o = P + (ps - P).regular() * c.r;
        return c;
    }
};

double disToSegment(Point s, Point a, Point b)  //点到线段ab最短距离
{
    if (a == b) return (s - a).length();
    Vector v1 = b - a, v2 = s - a, v3 = s - b;
    if (sgn(v1 * v2) < 0) return v2.length();
    if (sgn(v1 * v3) > 0) return v3.length();
    return Line(a, b - a).disToPoint(s);
}
double area(Point a, Point b, Point c)  //三角形面积
{
    return fabs((b - a).cross(c - a) / 2);
}
bool onSegment(Point p, Point a, Point b)  //点是否在线段上
{
    return (sgn((a - p).cross(b - p)) == 0 && sgn((a - p) * (b - p)) < 0) || p == a || p == b;
}
int insWithSegment(Line l, Point a, Point b, Point& ans)  //直线与线段相交
{
    if (l.parrelTo(Line(a, b - a))) return -1;  //不相交
    ans = l.insWithLine(Line(a, b - a));
    if (!onSegment(ans, a, b)) return -1;
    if (ans == a || ans == b) return 0;  //与端点相交
    return 1;                            //相交
}
bool segmentIns(Point a1, Point a2, Point b1, Point b2)  //判断线段是否相交
{
    double c1 = (a2 - a1).cross(b1 - a1), c2 = (a2 - a1).cross(b2 - a1);
    double c3 = (b2 - b1).cross(a1 - b1), c4 = (b2 - b1).cross(a2 - b1);
    if (!sgn(c1) || !sgn(c2) || !sgn(c3) || !sgn(c4)) {  //控制是否可以在顶点处相交
        bool f1 = onSegment(b1, a1, a2);
        bool f2 = onSegment(b2, a1, a2);
        bool f3 = onSegment(a1, b1, b2);
        bool f4 = onSegment(a2, b1, b2);
        bool f = (f1 | f2 | f3 | f4);
        return f;
    }
    return sgn(c1) * sgn(c2) < 0 && sgn(c3) * sgn(c4) < 0;
}
double polygonArea(Point* p, int n)  //计算多边形有向面积,点需要按顺序排列
{
    double s = 0;
    for (int i = 1; i < n - 1; i++) s += (p[i] - p[0]).cross(p[i + 1] - p[0]) / 2;
    return s;
}
double polygonPerimeter(Point* p, int n)  //计算多边形周长,点需要按顺序排列
{
    double ans = 0;
    p[n] = p[0];
    for (int i = 0; i < n; i++) ans += p[i].dist(p[i + 1]);
    return ans;
}
int quad(Point a)  // 判断象限的函数，每个象限包括半个坐标轴
{
    if (cmp(a.x, 0) > 0 && cmp(a.y, 0) >= 0) return 1;
    if (cmp(a.x, 0) <= 0 && cmp(a.y, 0) > 0) return 2;
    if (cmp(a.x, 0) < 0 && cmp(a.y, 0) <= 0) return 3;
    if (cmp(a.x, 0) >= 0 && cmp(a.y, 0) < 0) return 4;
}

void sortByPolarAngle(Point at, Point* begin, Point* end)  //极角逆时针排序，以at为极点
{
    sort(begin, end, [&](Point a, Point b) {
        a = a - at, b = b - at;
        int l1 = quad(a), l2 = quad(b);
        if (l1 == l2) {
            double c = a.cross(b);
            return cmp(c, 0) > 0 || (cmp(c, 0) == 0 && a.length() < b.length());
        }
        return l1 < l2;
    });
}

int gramhamScan(Point* res, int n, Point* ans)  //求凸包,O(nlogn),ans存放答案,返回凸包上点的数量
{
    if (n < 1) return 0;
    Point p = res[0];
    int k = 0, top = 1;
    for (int i = 1; i < n; i++) {
        if (res[i] < p) p = res[i], k = i;
    }
    swap(res[0], res[k]), sortByPolarAngle(res[0], res + 1, res + n), ans[0] = res[0];
    for (int i = 1; i < n; i++) {
        while (top > 2 && sgn((ans[top - 1] - ans[top - 2]).cross(res[i] - ans[top - 2])) <= 0) --top;
        ans[top++] = res[i];
    }
    return top;
}
double rotatingCalipers(Point* poly, int n)  //旋转卡壳,返回凸包直径
{
    double ans = 0;
    int j = 2;
    poly[n] = poly[0];
    for (int i = 0; i < n; i++) {
        while (area(poly[i], poly[i + 1], poly[j]) < area(poly[i], poly[i + 1], poly[j + 1])) {
            ++j;
            if (j > n) j = 0;
        }
        ans = max(ans, max(poly[i].dist(poly[j]), poly[i + 1].dist(poly[j])));
    }
    return ans;
}
int isPointInPolygon(Point p, Point* poly, int n)  //点在多边形内1;在外-1;在边界上0
{
    int wn = 0;
    for (int i = 0; i < n; ++i) {
        if (onSegment(p, poly[i], poly[(i + 1) % n])) return 0;
        int k = sgn((poly[(i + 1) % n] - poly[i]).cross(p - poly[i]));
        int d1 = sgn(poly[i].y - p.y);
        int d2 = sgn(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 -1;
}