计算几何模板

学校大佬写的，本人对某些错误做修改后发发出
如果后续使用过程发现错误或需要补充的地方，欢迎下方评论区指正，我会第一时间修改

namespace Computational_Geomrtry //计算几何命称空间
{
#include <algorithm>
#include <cmath>
#include <vector>
#define PI acos(-1)
    const double eps(1e-10);                                      //精度
    inline bool is0(double x) { return (x > 0 ? x : -x) <= eps; } //判断是否为0
    inline bool equal(double x, double y) { return is0(x - y); }  //判断是否相等

    class Point //点
    {
    public:
        Point() {}
        Point(double _x, double _y) : x(_x), y(_y) {}

    public:
        double x, y;

    public:
        bool operator==(const Point &b) { return is0(x - b.x) && is0(y - b.y); }
        void operator=(const Point &b) { x = b.x, y = b.y; }
    };
    const Point O(0.0, 0.0);    //坐标原点
    const Point VecX(1.0, 0.0); //x轴单位向量
    const Point VecY(0.0, 1.0); //y轴单位向量
    //点a到点b距离
    inline double len(const Point a, const Point b)
    {
        double dx(a.x - b.x), dy(a.y - b.y);
        return sqrt(dx * dx + dy * dy);
    }

    class Vector //向量
    {
    public:
        Vector() {}
        Vector(Point a) : x(a.x), y(a.y) {}
        Vector(Point a, Point b) : x(b.x - a.x), y(b.y - a.y) {}
        Vector(double _x, double _y) : x(_x), y(_y){};

    public:
        double x, y;

    public:
        double operator*(const Vector &b) { return x * b.x + y * b.y; }        //点积
        double operator^(const Vector &b) { return x * b.y - b.x * y; }        //叉积:右手定则判断正负
        Vector operator-(const Vector &b) { return Vector(x - b.x, y - b.y); } //向量相减
        Vector operator+(const Vector &b) { return Vector(x + b.x, y + b.y); } //向量相加

    public:
        double mold() const { return sqrt(x * x + y * y); }                                                         //向量模长
        double degree() const { return atan2(y, x); }                                                               //向量与x轴正半轴的夹角
        Vector Rotate(double seta) { return Vector{x * cos(seta) - y * sin(seta), x * sin(seta) + y * cos(seta)}; } ///逆时针旋转seta弧度,顺时针取负即可
    };

    class Circle //圆
    {
    public:
        Circle() {}
        Circle(const double _r) : o(O), r(_r) {}
        Circle(const Point &_o) : o(_o), r(0) {}
        Circle(const Point &_o, const double _r) : o(_o), r(_r) {}
        Circle(const Point &_o, const Point &_p) : o(_o), r(len(_o, _p)) {}

    public:
        Point o;  //圆心
        double r; //半径
    };

    /*****************************************************/

    //点p到线段<a,b>距离
    inline double pointToSegment(const Point &a, const Point &b, const Point &p)
    {
        double ap_ab(Vector(a, b) * Vector(a, p));
        if (ap_ab <= 0.0)
            return len(a, p);
        double d2(len(a, b));
        d2 *= d2;
        if (ap_ab >= d2)
            return len(b, p);
        double r(ap_ab / d2);
        double px(a.x + (b.x - a.x) * r);
        double py(a.y + (b.y - a.y) * r);
        return len(p, Point(px, py));
    }

    //点p到直线<a,b>距离
    inline double pointToLine(const Point &a, const Point &b, const Point &p)
    {
        double ap_ab(fabs(Vector(a, p) * Vector(a, b)));
        double d(ap_ab / len(a, b)), l(len(a, p));
        return sqrt(l * l - d * d);
    }

    //求<a,b,c>三点组成三角形的面积，若返回结果为0，则三点共线
    inline double getS(const Point &a, const Point &b, const Point &c)
    {
        return abs(0.5 * (Vector(a, b) ^ Vector(a, c)));
    }

    //判断线段[注意是线段不是直线!!!]<a,b>和<c,d>是否相交
    inline bool isIntersect(const Point &a, const Point &b, const Point &c, const Point &d)
    {
        Vector ac(a, c), ad(a, d), bc(b, c), bd(b, d);
        Vector ca(c, a), da(d, a), cb(c, b), db(d, b);
        return (ac ^ ad) * (bc ^ bd) <= 0.0 && (ca ^ cb) * (da ^ db) <= 0.0;
    }

    //判断两向量(直线)[<a,b>,<c,d>]是否平行,共线的话,返回-1,平行返回1,相交返回0
    inline int isParallel(const Point &a, const Point &b, const Point &c, const Point &d)
    {
        if (is0(getS(a, b, c)) && is0(getS(a, b, d)))
            return -1;
        return is0(Vector(a, b) ^ Vector(c, d));
    }

    //判断两线段(直线,向量)[<a,b>,<c,d>]是否垂直:点积是不是0
    inline bool isVertical(const Point &a, const Point &b, const Point &c, const Point &d)
    {
        return is0((Vector(a, b) * Vector(c, d)));
    }

    /*
直线的一般方程为F(x)=ax+by+c=0.既然我们已经知道直线的两个点,假设为(x0,y0),(x1,y1)那么可以得到a=y0–y1,b=x1–x0,c=x0y1–x1y0.
    */

    //求出<a,b> <c,d> 的交点坐标
    inline Point getCrossPoint(const Point &a, const Point &b, const Point &c, const Point &d)
    {
        double a0(a.y - b.y), b0(b.x - a.x), c0(a.x * b.y - b.x * a.y);
        double a1(c.y - d.y), b1(d.x - c.x), c1(c.x * d.y - d.x * c.y);
        double D(a0 * b1 - a1 * b0);
        if (is0(D))
            return Point(1ull << 63, 1ull << 63);
        return Point((b0 * c1 - b1 * c0) / D, (a1 * c0 - a0 * c1) / D);
    }

    //求三角形<a,b,c>的外心(外接圆的圆心)[外接圆的半径:a*b*c*0.25/getS(a,b,c)]
    inline Circle getCircumcircle(const Point &a, const Point &b, const Point &c)
    {
        double A1(2 * (b.x - a.x)), B1(2 * (b.y - a.y)), C1(b.x * b.x - a.x * a.x + b.y * b.y - a.y * a.y);
        double A2(2 * (c.x - b.x)), B2(2 * (c.y - b.y)), C2(c.x * c.x - b.x * b.x + c.y * c.y - b.y * b.y);
        double D(A1 * B2 - A2 * B1);
        Point ans((C1 * B2 - C2 * B1) / D, (A1 * C2 - A2 * C1) / D);
        return Circle(ans, len(ans, a));
    }

    //求三角形<a,b,c>内接圆[内切圆半径:2*getS(a,b,c)/(len(a,b)+len(b,c)+len(a,c))]
    inline Circle InscribedCircle(const Point &a, const Point &b, const Point &c)
    {
        double A(len(b, c)), B(len(a, c)), C(len(a, b)), sum(A + B + C);
        Point ans((A * a.x + B * b.x + C * c.x) / sum, (A * a.y + B * b.y + C * c.y) / sum);
        return Circle(ans, 2.0 * getS(a, b, c) / (A + B + C));
    }

    //最小圆覆盖(Minimum circle cover):求出到包含所有点的最小圆[精度有些问题]
    inline Circle MCC(Point *pbegin, int n)
    {
        Point *&p(pbegin);
        std::random_shuffle(p, p + n);
        Circle ans;
        for (int i(0); i < n; ++i)
        {
            if (len(ans.o, p[i]) > ans.r)
            {
                ans = Circle(p[i]);
                for (int j(0); j < i; ++j)
                {
                    if (len(ans.o, p[j]) > ans.r)
                    {
                        ans.o.x = (p[i].x + p[j].x) * 0.5;
                        ans.o.y = (p[i].y + p[j].y) * 0.5;
                        ans.r = len(p[i], p[j]) * 0.5;
                        for (int k(0); k < j; ++k)
                            if (len(ans.o, p[k]) > ans.r)
                                ans = getCircumcircle(p[i], p[j], p[k]);
                    }
                }
            }
        }
        return ans;
    }

    // 极角排序，以x轴正半轴作为基准
    void PA_Sort(Vector *pbegin, Vector *pend)
    {
        auto cmp = [&](const Vector &a, const Vector &b)
        {
            if (!is0(atan2(a.y, a.x) - atan2(b.y, b.x)))
                return atan2(a.y, a.x) < atan2(b.y, b.x);
            return a.mold() < b.mold();
        };
        std::sort(pbegin, pend, cmp);
    }
};