[学习笔记]从零开始的计算几何

\[\newcommand\fab[1]{{\left|{#1}\right|}} \newcommand\vec[1]{{\overrightarrow{#1}}} \]

目录

〇、前言 ¶

写算几的题的时候，注意精度问题，不然会精尽人亡。

壹、基础功能 ¶

由于 \(\tt double\) 的存储本身就有丢精问题，所以我们的程序中，所有涉及 \(\tt double\) 的判断应该都有一个误差容许值 \(\epsilon\)，并认为若 \(\fab{a-b}\le \epsilon\) 则判定 \(a=b\)，在我所有的代码中，如果没有特殊要求，一般都是 \(\epsilon=10^{-8}\).

然后就有以下比较器：

/** <---------------------- Basic function ----------------------> */
const double eps=1e-8;
const double Pi=acos(-1.0);
inline int comper(const double x, const double y){
    if(fab(x-y)<eps) return 0;
    return x>y? 1: -1;
}

贰、点 & 向量 ¶

◆ 向量

先说向量，向量需要实现的功能肯定包括加减，内外积，模长以及旋转，于是就有了这个板子：

struct Vec2d{
    double x, y;
    Vec2d(){}
    Vec2d(double X, double Y): x(X), y(Y){}
    inline Vec2d operator -(const Vec2d rhs) const{ return Vec2d(x-rhs.x, y-rhs.y); }
    inline Vec2d operator +(const Vec2d rhs) const{ return Vec2d(x+rhs.x, y+rhs.y); }
    // dot multiplication
    inline double operator ^(const Vec2d rhs) const{ return x*rhs.x+y*rhs.y; }
    // cross multiplication
    inline double operator *(const Vec2d rhs) const{ return x*rhs.y-y*rhs.x; }
    inline Vec2d operator *(const double k) const{ return Vec2d(k*x, k*y); }
    inline Vec2d rotate(const double theta) const{
        return Vec2d(x*cos(theta)-y*sin(theta), x*sin(theta)+y*cos(theta));
    }
    inline double norm() const{ return sqrt((*this)^(*this)); }
    inline friend double norm(const Vec2d v){ return v.norm(); }
    inline void print() const{ printf("(%.5f, %.5f)\n", x, y); }
    inline friend void print(const Vec2d v){ v.print(); }
};
/** get the vertical vector */
inline Vec2d getVertical(const Vec2d v){ return Vec2d(v.y, -v.x); }

此外，还实现了一个得到垂直向量的函数。

对于二维向量外积的特殊说明：

对于向量 \(\vec a,\vec b\)，有 \(\fab{\vec a\times \vec b}\) 表示 \(a,b\) 向量围成的四边形面积。

如果没有绝对值：若 \(\vec a\times \vec b<0\)，则说明 \(\vec a\) 在 \(\vec b\) 的逆时针方向，反之为顺时针，若 \(\vec a\times \vec b=0\)，则 \(\vec a,\vec b\) 线性相关。

对于三维向量外积的特殊说明：

对于 \(\vec a,\vec b,\vec c\) 三个三维向量，若 \(\vec c=\vec a\times \vec b\).

记三个矩阵：

\[D_x= \begin{bmatrix} a.y&a.z \\ b.y&b.z \end{bmatrix} \quad D_y= \begin{bmatrix} a.z&a.x \\ b.z&b.x \end{bmatrix} \quad D_z= \begin{bmatrix} a.x&a.y \\ b.x&b.y \end{bmatrix} \]

那么 \(\vec c=(\det D_x,\det D_y,\det D_z)\)，且 \(\fab{\vec c}=\fab{\vec a}\fab{\vec b}\sin\lang \vec a,\vec b\rang\).

◆ 点

点就比向量好很多了，注意两个点减出来会是一个向量。

struct Point{
    double x, y;
    Point(){}
    Point(double X, double Y): x(X), y(Y){}
    inline Point operator +(const Vec2d v) const{ return Point(x+v.x, y+v.y); }
    inline Vec2d operator -(const Point rhs) const{ return Vec2d(x-rhs.x, y-rhs.y); }
    inline void print() const{ printf("(%.5f, %.5f)\n", x, y); }
    inline Vec2d ptov() const{ return Vec2d(x, y); }
};
inline Vec2d ptov(const Point p){ return Vec2d(p.x, p.y); }
inline Point vtop(const Vec2d v){ return Point(v.x, v.y); }
inline bool same(const Point p1, const Point p2){
    return comper(p1.x, p2.x)==0 && comper(p1.y, p2.y)==0;
}

其中还实现了两个向量与点互转的函数以及两个点相同的比较。

贰、直线 & 线段 ¶

◆ 直线

先说直线，线段比直线更复杂。

直线的定义

我们采用直线的参数表示法，即找一个点，再找一个方向向量。

struct Line{
    Point p; Vec2d v;
    Line(){}
    Line(Point P, Vec2d V): p(P), v(V){}
    inline void print(){
        printf("point : "); p.print();
        printf("vecto : "); v.print();
    }
};

点到直线的距离

直接用平行四边形面积除以底边长，就可以得到高了。

inline double getDis(const Line l, const Point p){
    Vec2d v=p-l.p; double siz=fab(v*l.v);
    return siz/l.v.norm();
}

点是否在直线上

判断距离与 \(0\) 即可。

inline bool isOnLine(const Line l, const Point p){
    double dis=getDis(l, p);
    return comper(dis, 0)==0;
}

判断两直线斜率是否相同

判断方向向量是否线性相关即可。

inline bool parallel(Line a, Line b){
    double c=a.v*b.v;
    if(comper(c, 0)==0) return 1;
    return 0;
}

求两直线的交点

这是个大话题，我们不妨先解个方程，先设

\[l_1:p=p_1+t\vec{v_1}\quad l_2:p=p_3+s\vec{v_2} \]

其中 \(\vec{v_1}=p_2-p_1,\vec{v_2}=p_4-p_3\)

将 \(x,y\) 分离，得到

\[t(x_2-x_1)-s(x_4-x_3)=x_3-x_1 \\ t(y_2-y_1)-s(y_4-y_3)=y_3-y_1 \]

使用克莱姆法则求解该方程，于是，有

\[D_0= \fab{ \begin{matrix} x_2-x_1&-(x_4-x_3) \\ y_2-y_1&-(y_4-y_3) \end{matrix} } = -\fab{ \begin{matrix} x_2-x_1&x_4-x_3 \\ y_2-y_1&y_4-y_3 \end{matrix} } \\ D_t= \fab{ \begin{matrix} x_3-x_1&-(x_4-x_3) \\ y_3-y_1&-(y_4-y_3) \end{matrix} } = \fab{ \begin{matrix} x_1-x_3&x_4-x_3 \\ y_1-y_3&y_4-y_3 \end{matrix} } \\ D_s= \fab{ \begin{matrix} x_2-x_1&x_3-x_1 \\ y_2-y_1&y_3-y_1 \end{matrix} } = -\fab{ \begin{matrix} x_2-x_1&x_1-x_3 \\ y_2-y_1&y_1-y_3 \end{matrix} } \]

经过整理之后发现这有点像叉乘，于是我们可以定义

\[\vec a=p_2-p_1=\vec{v_1}\quad \vec b=p_4-p_3=\vec{v_2}\quad \vec v=p_1-p_3 \]

那么就有

\[D_0=-\vec{v_1}\times\vec{v_2}\quad D_t=\vec v\times \vec{v_2}\quad D_s=-\vec{v_1}\times v \]

所以，有

\[t={D_t\over D_0}={\vec{v_2}\times \vec v\over \vec{v_1}\times \vec{v_2}} \\ s={D_s\over D_0}={\vec{v_1}\times \vec v\over \vec{v_1}\times \vec{v_2}} \]

这个结构很对称诶。解出来之后代入直线方程就行了。

/** @warning @p a and @p b shouldn't be the same Line or parallel */
inline Point getIntersect(Line a, Line b){
    double d0=a.v*b.v;
    assert(comper(d0, 0)!=0);
    double t=b.v*(a.p-b.p)/d0;
    return a.p+a.v*t;
}

◆ 线段

线段的定义

存下两个端点，这应该是最方便的。另外还实现了一个得到该线段所在直线的函数。

struct Segmt{
    Point p0, p1;
    Segmt(){}
    Segmt(const Point P0, const Point P1): p0(P0), p1(P1){}
    inline Point& operator [](const int i){
        assert(i==0 || i==1); return i==0? p0: p1;
    }
    inline Line SegtoLine(){ return Line(p0, p1-p0); }
    inline void print(){
        p0.print(); p1.print();
    }
};

线段相交

跨立实验：判断 \(p_1,p_2\) 是否在直线 \(p_3p_4\) 的两侧。再判断 \(p_3,p_4\) 是否在直线 \(p_1p_2\) 的两侧。若通过跨立实验，则说明线段有交点。

但还应需要考虑两线段共线的情况，以及一条线段的端点在另一条线段上的情况。

inline bool segSegIntersect(Segmt a, Segmt b){
    for(int i=0; i<2; ++i) if(isOnSeg(b, a[i])) return 1;
    for(int j=0; j<2; ++j) if(isOnSeg(a, b[j])) return 1;
    if(parallel(a, b)) return 0;
    double c0=(a[1]-a[0])*(b[0]-a[0]), c1=(a[1]-a[0])*(b[1]-a[0]);
    double c2=(b[1]-b[0])*(a[0]-b[0]), c3=(b[1]-b[0])*(a[1]-b[0]);
    if(comper(c0, 0)*comper(c1, 0)>0 || comper(c2, 0)*comper(c3, 0)>0)
        return 0;
    return 1;
}

点到线段的距离

先判断垂足是否落在线段上即可。

inline double getDis(Segmt l, const Point p){
    double d0=(p-l[0])^(l[1]-l[0]), d1=(p-l[1])^(l[0]-l[1]);
    if(comper(d0, 0)>=0 && comper(d1, 0)>=0)
        return getDis(Line(l[0], l[1]-l[0]), p);
    return comper(d0, 0)<0? norm(p-l[0]): norm(p-l[1]);
}

点是否在线段上

判断距离。

inline bool isOnSeg(const Segmt l, const Point p){
    return comper(getDis(l, p), 0)==0;
}

线段是否与直线交

缩略版跨立实验。

inline bool segLineIntersect(Segmt s, const Line l){
    if(isOnLine(l, s[0]) || isOnLine(l, s[1])) return 1;
    double c0=l.v*(s[0]-l.p), c1=l.v*(s[1]-l.p);
    if(comper(c0, 0)*comper(c1, 0)<0) return 1;
    return 0;
}

线段斜率比较

比所在直线的就行了。

inline bool parallel(Segmt a, Segmt b){
    double c=(a[1]-a[0])*(b[1]-b[0]);
    if(comper(c, 0)==0) return 1;
    return 0;
}

线段的中垂线

中点+垂直向量即可。

inline Line getVerticalLine(Segmt s){
    return Line(getMid(s[0], s[1]), getVertical(s[0]-s[1]));
}

叁、多边形 & 凸包 ¶

◆ 多边形的表示

注意顶点顺序，拿一个 \(\tt vector\) 来存就可以了。

多边形面积

对多边形进行三角剖分，用叉乘结果除以 \(2\) 就可以得到面积了。注意，应当直接使用有向面积。该方法适用于复杂多边形。

inline double getArea(vector<Point>poly){
    int n=poly.size(); double s=0;
    for(int i=0; i<n; ++i) s+=poly[i].ptov()*poly[(i+1)%n].ptov();
    return s*0.5;
}

判断一个点是否在多边形中

由于另一种方法不适用与复杂多边形，故而此处不作介绍。此处使用射线法判断：从该点引出一条射线，看射线与多边形交点数量，若为奇数，则在内部中，否则不在。

// whether vector @p x is in the middle of @p a and @p b
inline bool isInMiddle(Vec2d a, Vec2d b, Vec2d x){
    if(comper(a*b, 0)<0) swap(a, b);
    return comper(a*x, 0)>=0 && comper(b*x, 0)<=0;
}
inline bool isInside(vector<Point>poly, Point p){
    int n=poly.size();
    rep(_, 1, 5){
        double angle=(rand()/(double(RAND_MAX))-0.5)*2*Pi;
        Vec2d v(cos(angle), sin(angle));
        bool res=0;
        for(int i=0; i<n; ++i){
            if(isOnSeg(Segmt(poly[i], poly[(i+1)%n]), p))
                return 1;
            if(isInMiddle(poly[i]-p, poly[(i+1)%n]-p, v))
                res=!res;
        }
        if(res) return 1;
    }
    return 0;
}

找到重心

先对多边形进行三角剖分，得到每个三角的重心，然后再代权平均即可。

inline Point getCG(vector<Point>poly){
    double area=0;
    int n=poly.size();
    Point ret(0, 0);
    for(int i=0; i<n; ++i){
        double t=poly[i].ptov()*poly[(i+1)%n].ptov();
        area+=t;
        ret.x+=(poly[i].x+poly[(i+1)%n].x)*t;
        ret.y+=(poly[i].y+poly[(i+1)%n].y)*t;
    }
    ret.x/=3, ret.x/=area;
    ret.y/=3, ret.y/=area;
    return ret;
}

◆ 凸包

定义就不说了，还是用一个 \(\tt vector\) 来存就可以了。

求多个点的凸包

\(\rm Graham\) 扫描法

找到最左下角的点，记为 \(p_0\)，并以其为原点，对剩下的点进行极角排序，极角相同的，距离小的在前。

先加入最开始的三个点（包含 \(p_0\)），然后按极角序枚举其他点，看新点加入时，是在上一个点基础上左拐，则加入栈，否则一直弹，直到左拐或者剩下 \(p_0\) 为止。

该算法主要问题在于，不能完整地保存凸包边上的所有点（但是应该可以再暴力扫一遍）。

Point sta[STACK_SIZE+5]; int ed;
inline vector<Point> Graham(vector<Point>poly){
    Point pivot; int k=0, n=poly.size();
    if(n==1) return {poly[0]};
    for(int i=1; i<n; ++i)
        if(comper(poly[i].y, poly[k].y)<0 || (comper(poly[i].y, poly[k].y)==0 && comper(poly[i].x, poly[k].x)<0))
            k=i;
    pivot=poly[k];
    swap(poly[0], poly[k]);
    auto cmp=[&](const Point a, const Point b){
        Vec2d v1=a-pivot, v2=b-pivot;
        double c=v1*v2;
        return comper(c, 0)>0 || (comper(c, 0)==0 && comper(norm(v1), norm(v2))<0);
    };
    sort(poly.begin()+1, poly.end(), cmp);
    sta[0]=poly[0], sta[1]=poly[1], ed=1;
    for(int i=2; i<n; ++i){
        while(ed && comper((sta[ed]-sta[ed-1])*(poly[i]-sta[ed]), 0)<=0)
            --ed;
        sta[++ed]=poly[i];
    }
    vector<Point>convex;
    for(int i=0; i<=ed; ++i)
        convex.push_back(sta[i]);
    return convex;
}

对 \(\rm Graham\) 扫描法的极角排序加以改进：

把按照坐标关键字排序后的前两个点 \(p_0,p_1\) 入栈，对于后面的每个点 \(p_2,\cdots,p_{n-1}\)

若栈中有两个元素及以上，设栈顶为 \(k_2\)，栈顶下一个为 \(k_1\)，考虑 \(k_1-k_2-p_i\) 是否向左拐，如果向右拐或者直走，就弹出 \(k_2\) 直到不再弹栈或者栈中只有一个元素时，把 \(p_i\) 入栈。

这样得到点集的下凸壳，再把排序的点倒过来，做一遍得到点集的上凸壳。上下拼在一起，得到完整的凸包。

Minkowski Sum

在两个按照极角从小到大有序的凸包中进行合并。

假设现在做到 \(A[i]\) 和 \(B[j]\) 两个点，把 \(A[i]+B[j]\) 放入集合，然后用叉积比较 \(A[i+1]-A[i],B[j+1]-B[j]\) 的极角序，把极角序较小者的指针 \(+1\)，如果极角序相同则都 \(+1\).

inline vector<Point> Minkowski_sum(vector<Point>p1, vector<Point>p2){
    vector<Point>ret;
    int n1=p1.size(), n2=p2.size();
    int i=0, j=0;
    while(i<n1 && j<n2){
        Vec2d a=p1[(i+1)%n1]-p1[i], b=p2[(j+1)%n2]-p2[j];
        ret.push_back(vtop(p1[i].ptov()+p2[j].ptov()));
        if(comper(a*b, 0)==0) ++i, ++j;
        else if(comper(a*b, 0)<0) ++j;
        else ++i;
    }
    while(i<n1) ret.push_back(vtop(p1[i++].ptov()+p2[0].ptov()));
    while(j<n2) ret.push_back(vtop(p1[0].ptov()+p2[j++].ptov()));
    return ret;
}

肆、其他东西 ¶

◆ 圆

圆的存储

一个点，一个半径，就完事了。

struct Circle{
    Point c; double r;
    Circle(){}
    Circle(const Point C, const double R): c(C), r(R){}
    inline double area(const double para=1){ return Pi*r*r*para; }
    inline Point& operator [](const int i){
        assert(i==0); return c;
    }
};
inline double getArea(Circle C, const double para=1){ return C.area(para); }

在里面还实现了一个得到面积的函数。

三点定圆

垂线相交即可，注意特判三点共线。

inline Circle getCC(const Point a, const Point b, const Point c){
    if(comper((b-a)*(c-b), 0)==0){
        double dis=0; Point x, y;
        if(dis<norm(a-b)) dis=norm(a-b), x=a, y=b;
        if(dis<norm(a-c)) dis=norm(a-c), x=a, y=c;
        if(dis<norm(b-c)) dis=norm(b-c), x=b, y=c;
        return Circle(getMid(x, y), dis*0.5);
    }
    Line l1=getVerticalLine(Segmt(a, b)), l2=getVerticalLine(Segmt(c, b));
    Point p=getIntersect(l1, l2);
    return Circle(p, norm(a-p));
}

点圆关系

判断距离与半径就行了。

inline int isInside(Circle C, Point p){
    int res=comper(norm(C[0]-p), C.r);
    if(res<0) return 1;
    if(res==0) return 0;
    return -1;
}

最小圆覆盖

使用随机增量法。

先让 \(Di=\{p_1,p_i\}\)，逐个判断点 \(\{p_2,p_3,\cdots ,p_{i-1}\}\)，如果存在 \(p_j\notin D_i\)，则 \(D_i=\{p_j,p_i\}\).

同时，再依次判断点 \(\{p_1,p_2,\cdots ,p_{j-1}\}\)，如果存在 \(p_k\notin D_i\)，则 \(D_i=\{p_k,p_j,p_i\}\).

因为三点唯一地确定一个圆，所以只需在此基础上判断其他的点是否位于此圆内，不停地更新 \(p_k\).

当 \(k\) 的循环完成时，退出循环，更新 \(p_j\)；当 \(j\) 的循环结束时，退出循环，更新 \(p_i\). 当 \(i=n\) 时，算法完成。

将所有点随机化之后，由于一个圆由三个点确定，所以触发 \(i,j\) 循环的概率是 \(3\over i\)，同理，\(j,k\) 循环概率为 \(3\over j\)，所以总复杂度期望线性。

inline Circle getMinCirCover(vector<Point>p){
    random_shuffle(p.begin(), p.end());
    int n=p.size();
    Circle C(p[0], 0);
    for(int i=1; i<n; ++i){
        if(isInside(C, p[i])<0){
            C=Circle(p[i], 0);
            for(int j=0; j<i; ++j)
                if(isInside(C, p[j])<0){
                    C[0]=getMid(p[i], p[j]), C.r=norm(p[i]-C[0]);
                    for(int k=1; k<j; ++k) if(isInside(C, p[k])<0)
                        C=getCC(p[i], p[j], p[k]);
                }
        }
    }
    return C;
}

圆是否包含

判断圆心距与半径差即可。

inline bool isContain(Circle in, Circle out){
    if(comper(out.r-norm(out[0]-in[0]), in.r)>=0)
        return 1;
    return 0;
}

圆与圆相交的面积

大体思路是用两个扇形面积减去一个四边形面积。

先用余弦求出圆心角，再用正弦求出四边形面积，作差即可，注意判断包含与相离的情况。

inline double crossArea(Circle C1, Circle C2){
    if(isContain(C1, C2) || isContain(C2, C1))
        return min(C1.area(), C2.area());
    else if(comper(norm(C1[0]-C2[0]), C1.r+C2.r)>=0)
        return 0;
    double d=norm(C1[0]-C2[0]);
    double alpha=acos((C1.r*C1.r+d*d-C2.r*C2.r)/(2*C1.r*d));
    double beta=acos((C2.r*C2.r+d*d-C1.r*C1.r)/(2*C2.r*d));
    return C1.area(alpha/Pi)+C2.area(beta/Pi)-C1.r*d*sin(alpha);
}

伍、完整代码 ¶

// #define NDEBUG
#include<cassert>

namespace Elaina{
    #define rep(i, l, r) for(int i=(l), i##_end_=(r); i<=i##_end_; ++i)
    #define drep(i, l, r) for(int i=(l), i##_end_=(r); i>=i##_end_; --i)
    #define fi first
    #define se second
    #define mp(a, b) make_pair(a, b)
    #define Endl putchar('\n')
    #define mmset(a, b) memset(a, b, sizeof a)
    // #define int long long
    typedef long long ll;
    typedef unsigned long long ull;
    typedef pair<int, int> pii;
    typedef pair<ll, ll> pll;
    template<class T>inline T fab(T x){ return x<0? -x: x; }
    template<class T>inline void getmin(T& x, const T rhs){ x=min(x, rhs); }
    template<class T>inline void getmax(T& x, const T rhs){ x=max(x, rhs); }
    template<class T>inline T readin(T x){
        x=0; int f=0; char c;
        while((c=getchar())<'0' || '9'<c) if(c=='-') f=1;
        for(x=(c^48); '0'<=(c=getchar()) && c<='9'; x=(x<<1)+(x<<3)+(c^48));
        return f? -x: x;
    }
    template<class T>inline void writc(T x, char s='\n'){
        static int fwri_sta[1005], fwri_ed=0;
        if(x<0) putchar('-'), x=-x;
        do fwri_sta[++fwri_ed]=x%10, x/=10; while(x);
        while(putchar(fwri_sta[fwri_ed--]^48), fwri_ed);
        putchar(s);
    }
}
using namespace Elaina;

namespace Geometry{

    /** <---------------------- Basic function ----------------------> */
    const double eps=1e-8;
    const double Pi=acos(-1.0);
    const int INSIDE_CHECK_DATA=5;
    const int STACK_SIZE=1e5;
    inline int comper(const double x, const double y){
        if(fab(x-y)<eps) return 0;
        return x>y? 1: -1;
    }

    /** <---------------------- Vector & Point ----------------------> */
    struct Vec2d{
        double x, y;
        Vec2d(){}
        Vec2d(double X, double Y): x(X), y(Y){}
        inline Vec2d operator -(const Vec2d rhs) const{ return Vec2d(x-rhs.x, y-rhs.y); }
        inline Vec2d operator +(const Vec2d rhs) const{ return Vec2d(x+rhs.x, y+rhs.y); }
        // dot multiplication
        inline double operator ^(const Vec2d rhs) const{ return x*rhs.x+y*rhs.y; }
        // cross multiplication
        inline double operator *(const Vec2d rhs) const{ return x*rhs.y-y*rhs.x; }
        inline Vec2d operator *(const double k) const{ return Vec2d(k*x, k*y); }
        inline Vec2d rotate(const double theta) const{
            return Vec2d(x*cos(theta)-y*sin(theta), x*sin(theta)+y*cos(theta));
        }
        inline double norm() const{ return sqrt((*this)^(*this)); }
        inline friend double norm(const Vec2d v){ return v.norm(); }
        inline void print() const{ printf("(%.5f, %.5f)\n", x, y); }
        inline friend void print(const Vec2d v){ v.print(); }
    };
    struct Point{
        double x, y;
        Point(){}
        Point(double X, double Y): x(X), y(Y){}
        inline Point operator +(const Vec2d v) const{ return Point(x+v.x, y+v.y); }
        inline Vec2d operator -(const Point rhs) const{ return Vec2d(x-rhs.x, y-rhs.y); }
        inline void print() const{ printf("(%.5f, %.5f)\n", x, y); }
        inline Vec2d ptov() const{ return Vec2d(x, y); }
    };
    inline Point getMid(const Point a, const Point b){
        return Point((a.x+b.x)/2, (a.y+b.y)/2);
    }
    inline Vec2d ptov(const Point p){ return Vec2d(p.x, p.y); }
    /** get the vertical vector */
    inline Vec2d getVertical(const Vec2d v){ return Vec2d(v.y, -v.x); }
    inline Point vtop(const Vec2d v){ return Point(v.x, v.y); }
    inline bool same(const Point p1, const Point p2){
        return comper(p1.x, p2.x)==0 && comper(p1.y, p2.y)==0;
    }

    /** <---------------------- Line & Segment ----------------------> */
    struct Line{
        Point p; Vec2d v;
        Line(){}
        Line(Point P, Vec2d V): p(P), v(V){}
        inline void print(){
            printf("point : "); p.print();
            printf("vecto : "); v.print();
        }
    };
    struct Segmt{
        Point p0, p1;
        Segmt(){}
        Segmt(const Point P0, const Point P1): p0(P0), p1(P1){}
        inline Point& operator [](const int i){
            assert(i==0 || i==1); return i==0? p0: p1;
        }
        inline Line SegtoLine(){ return Line(p0, p1-p0); }
        inline void print(){
            p0.print(); p1.print();
        }
    };
    inline bool parallel(Segmt a, Segmt b){
        double c=(a[1]-a[0])*(b[1]-b[0]);
        if(comper(c, 0)==0) return 1;
        return 0;
    }
    inline bool parallel(Line a, Line b){
        double c=a.v*b.v;
        if(comper(c, 0)==0) return 1;
        return 0;
    }
    inline double getDis(const Line l, const Point p){
        Vec2d v=p-l.p; double siz=fab(v*l.v);
        return siz/l.v.norm();
    }
    inline double getDis(Segmt l, const Point p){
        double d0=(p-l[0])^(l[1]-l[0]), d1=(p-l[1])^(l[0]-l[1]);
        if(comper(d0, 0)>=0 && comper(d1, 0)>=0)
            return getDis(Line(l[0], l[1]-l[0]), p);
        return comper(d0, 0)<0? norm(p-l[0]): norm(p-l[1]);
    }
    inline bool isOnLine(const Line l, const Point p){
        double dis=getDis(l, p);
        return comper(dis, 0)==0;
    }
    inline bool isOnSeg(const Segmt l, const Point p){
        return comper(getDis(l, p), 0)==0;
    }
    inline bool segSegIntersect(Segmt a, Segmt b){
        for(int i=0; i<2; ++i) if(isOnSeg(b, a[i])) return 1;
        for(int j=0; j<2; ++j) if(isOnSeg(a, b[j])) return 1;
        if(parallel(a, b)) return 0;
        double c0=(a[1]-a[0])*(b[0]-a[0]), c1=(a[1]-a[0])*(b[1]-a[0]);
        double c2=(b[1]-b[0])*(a[0]-b[0]), c3=(b[1]-b[0])*(a[1]-b[0]);
        if(comper(c0, 0)*comper(c1, 0)>0 || comper(c2, 0)*comper(c3, 0)>0)
            return 0;
        return 1;
    }
    inline bool segLineIntersect(Segmt s, const Line l){
        if(isOnLine(l, s[0]) || isOnLine(l, s[1])) return 1;
        double c0=l.v*(s[0]-l.p), c1=l.v*(s[1]-l.p);
        if(comper(c0, 0)*comper(c1, 0)<0) return 1;
        return 0;
    }
    inline bool coincident(Line a, Line b){
        if(!parallel(a, b)) return 0;
        if(comper(getDis(b, a.p), 0)==0) return 1;
        return 0;
    }
    /** @warning @p a and @p b shouldn't be the same Line or parallel */
    inline Point getIntersect(Line a, Line b){
        double d0=a.v*b.v;
        assert(comper(d0, 0)!=0);
        double t=b.v*(a.p-b.p)/d0;
        return a.p+a.v*t;
    }
    inline Line getVerticalLine(Segmt s){
        return Line(getMid(s[0], s[1]), getVertical(s[0]-s[1]));
    }

    /** <---------------------- polygon & convex ----------------------> */
    inline double getArea(vector<Point>poly){
        int n=poly.size(); double s=0;
        for(int i=0; i<n; ++i) s+=poly[i].ptov()*poly[(i+1)%n].ptov();
        return s*0.5;
    }
    // whether vector @p x is in the middle of @p a and @p b
    inline bool isInMiddle(Vec2d a, Vec2d b, Vec2d x){
        if(comper(a*b, 0)<0) swap(a, b);
        return comper(a*x, 0)>=0 && comper(b*x, 0)<=0;
    }
    inline bool isInside(vector<Point>poly, Point p){
        int n=poly.size();
        rep(_, 1, INSIDE_CHECK_DATA){
            double angle=(rand()/(double(RAND_MAX))-0.5)*2*Pi;
            Vec2d v(cos(angle), sin(angle));
            bool res=0;
            for(int i=0; i<n; ++i){
                if(isOnSeg(Segmt(poly[i], poly[(i+1)%n]), p))
                    return 1;
                if(isInMiddle(poly[i]-p, poly[(i+1)%n]-p, v))
                    res=!res;
            }
            if(res) return 1;
        }
        return 0;
    }
    Point sta[STACK_SIZE+5]; int ed;
    inline vector<Point> Graham(vector<Point>poly){
        Point pivot; int k=0, n=poly.size();
        if(n==1) return {poly[0]};
        for(int i=1; i<n; ++i)
            if(comper(poly[i].y, poly[k].y)<0 || (comper(poly[i].y, poly[k].y)==0 && comper(poly[i].x, poly[k].x)<0))
                k=i;
        pivot=poly[k];
        swap(poly[0], poly[k]);
        auto cmp=[&](const Point a, const Point b){
            Vec2d v1=a-pivot, v2=b-pivot;
            double c=v1*v2;
            return comper(c, 0)>0 || (comper(c, 0)==0 && comper(norm(v1), norm(v2))<0);
        };
        sort(poly.begin()+1, poly.end(), cmp);
        sta[0]=poly[0], sta[1]=poly[1], ed=1;
        for(int i=2; i<n; ++i){
            while(ed && comper((sta[ed]-sta[ed-1])*(poly[i]-sta[ed]), 0)<=0)
                --ed;
            sta[++ed]=poly[i];
        }
        vector<Point>convex;
        for(int i=0; i<=ed; ++i)
            convex.push_back(sta[i]);
        return convex;
    }
    inline Point getCG(vector<Point>poly){
        double area=0;
        int n=poly.size();
        Point ret(0, 0);
        for(int i=0; i<n; ++i){
            double t=poly[i].ptov()*poly[(i+1)%n].ptov();
            area+=t;
            ret.x+=(poly[i].x+poly[(i+1)%n].x)*t;
            ret.y+=(poly[i].y+poly[(i+1)%n].y)*t;
        }
        ret.x/=3, ret.x/=area;
        ret.y/=3, ret.y/=area;
        return ret;
    }
    inline vector<Point> Minkowski_sum(vector<Point>p1, vector<Point>p2){
        vector<Point>ret;
        int n1=p1.size(), n2=p2.size();
        int i=0, j=0;
        while(i<n1 && j<n2){
            Vec2d a=p1[(i+1)%n1]-p1[i], b=p2[(j+1)%n2]-p2[j];
            ret.push_back(vtop(p1[i].ptov()+p2[j].ptov()));
            if(comper(a*b, 0)==0) ++i, ++j;
            else if(comper(a*b, 0)<0) ++j;
            else ++i;
        }
        while(i<n1) ret.push_back(vtop(p1[i++].ptov()+p2[0].ptov()));
        while(j<n2) ret.push_back(vtop(p1[0].ptov()+p2[j++].ptov()));
        return ret;
    }
    
    /** <---------------------- others ----------------------> */
    struct Circle{
        Point c; double r;
        Circle(){}
        Circle(const Point C, const double R): c(C), r(R){}
        inline double area(const double para=1){ return Pi*r*r*para; }
        inline Point& operator [](const int i){
            assert(i==0); return c;
        }
    };
    inline double getArea(Circle C, const double para=1){ return C.area(para); }
    inline Circle getCC(const Point a, const Point b, const Point c){
        if(comper((b-a)*(c-b), 0)==0){
            double dis=0; Point x, y;
            if(dis<norm(a-b)) dis=norm(a-b), x=a, y=b;
            if(dis<norm(a-c)) dis=norm(a-c), x=a, y=c;
            if(dis<norm(b-c)) dis=norm(b-c), x=b, y=c;
            return Circle(getMid(x, y), dis*0.5);
        }
        Line l1=getVerticalLine(Segmt(a, b)), l2=getVerticalLine(Segmt(c, b));
        Point p=getIntersect(l1, l2);
        return Circle(p, norm(a-p));
    }
    /** @brief if the node is on the edge, return 0*/
    inline int isInside(Circle C, Point p){
        int res=comper(norm(C[0]-p), C.r);
        if(res<0) return 1;
        if(res==0) return 0;
        return -1;
    }
    inline Circle getMinCirCover(vector<Point>p){
        random_shuffle(p.begin(), p.end());
        int n=p.size();
        Circle C(p[0], 0);
        for(int i=1; i<n; ++i){
            if(isInside(C, p[i])<0){
                C=Circle(p[i], 0);
                for(int j=0; j<i; ++j)
                    if(isInside(C, p[j])<0){
                        C[0]=getMid(p[i], p[j]), C.r=norm(p[i]-C[0]);
                        for(int k=1; k<j; ++k) if(isInside(C, p[k])<0)
                            C=getCC(p[i], p[j], p[k]);
                    }
            }
        }
        return C;
    }
    inline bool isContain(Circle in, Circle out){
        if(comper(out.r-norm(out[0]-in[0]), in.r)>=0)
            return 1;
        return 0;
    }
    inline double crossArea(Circle C1, Circle C2){
        if(isContain(C1, C2) || isContain(C2, C1))
            return min(C1.area(), C2.area());
        else if(comper(norm(C1[0]-C2[0]), C1.r+C2.r)>=0)
            return 0;
        double d=norm(C1[0]-C2[0]);
        double alpha=acos((C1.r*C1.r+d*d-C2.r*C2.r)/(2*C1.r*d));
        double beta=acos((C2.r*C2.r+d*d-C1.r*C1.r)/(2*C2.r*d));
        return C1.area(alpha/Pi)+C2.area(beta/Pi)-C1.r*d*sin(alpha);
    }
}
using namespace Geometry;

陆、特别提醒 ¶

如果 \(\color{yellow}{\text{TLE}}\) 了，可以将 INSIDE_CHECK_DATA 调小，一般 \(1\) 也足够了。

如果 \(\color{purple}{\text{RE}}\) 了，看看是不是 STACK_SIZE 小了，或者是 assert() 被触发了。

如果 \(\color{red}{\text{WA}}\) 了，看看是不是精尽人亡了，可以将 double 换成 long double，或者交 \(\rm C++\) 而不是 \(\rm G++\).

posted @ 2021-08-06 22:42  Arextre  阅读(269)  评论(0)  编辑  收藏  举报