本文档整理自算法竞赛常用模板，包含计算几何领域的核心算法和数据结构。

基于叉积的极角排序

叉积定义

$$\text{cross}(\vec{a}, \vec{b}) = x_1 y_2 - y_1 x_2$$

几何意义

  • $> 0$：$\vec{b}$ 在 $\vec{a}$ 的逆时针方向
  • $< 0$：$\vec{b}$ 在 $\vec{a}$ 的顺时针方向
  • $= 0$：两向量共线

象限划分

象限 条件 顺序
0 $x > 0, y \ge 0$ 第一象限（含正y轴）
1 $x \le 0, y > 0$ 第二象限（含负x轴）
2 $x < 0, y \le 0$ 第三象限（含负y轴）
3 $x \ge 0, y < 0$ 第四象限（含正x轴）

极角排序模板

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
struct Point {
    long long x, y;
    Point(long long x = 0, long long y = 0) : x(x), y(y) {}
};

long long cross(const Point& a, const Point& b) {
    return a.x * b.y - a.y * b.x;
}

int quad(const Point& p) {
    if (p.x > 0 && p.y >= 0) return 0;
    if (p.x <= 0 && p.y > 0) return 1;
    if (p.x < 0 && p.y <= 0) return 2;
    return 3;
}

bool cmpPolar(const Point& a, const Point& b) {
    int qa = quad(a), qb = quad(b);
    if (qa != qb) return qa < qb;

    long long c = cross(a, b);
    if (c == 0) {  // 共线：按模长排序
        return a.x * a.x + a.y * a.y < b.x * b.x + b.y * b.y;
    }
    return c > 0;  // 逆时针方向极角小
}

相对于任意参考点的排序

1
2
3
4
5
6
Point ref;
vector<Point> trans;
for (auto& p : pts) {
    trans.push_back({p.x - ref.x, p.y - ref.y});
}
sort(trans.begin(), trans.end(), cmpPolar);

完整计算几何板子

以下是一套完整的计算几何板子，使用 long double 保证精度，涵盖点、线、线段、多边形、凸包、圆等核心几何结构及常用算法。

核心定义与工具函数

1
2
3
4
5
6
7
8
9
10
11
12
#include <bits/stdc++.h>
using namespace std;

using ld = long double;
constexpr ld EPS = 1e-8;
constexpr ld PI = 3.1415926535897932384l;
constexpr long long INF = 0x3f3f3f3f3f3f3f3f;

int sgn(ld x) {
    if (fabs(x) < EPS) return 0;
    return x < 0 ? -1 : 1;
}

点结构体

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
struct Point {
    ld x, y;
    Point() {}
    Point(ld xx, ld yy) : x(xx), y(yy) {}
    Point operator + (const Point& b) const { return Point(x + b.x, y + b.y); }
    Point operator - (const Point& b) const { return Point(x - b.x, y - b.y); }
    Point operator * (ld b) const { return Point(b * x, b * y); }
    Point operator / (ld b) const { return Point(x / b, y / b); }
    ld operator * (const Point& b) const { return x * b.x + y * b.y; }      // 点积
    ld operator % (const Point& b) const { return x * b.y - y * b.x; }      // 叉积
    bool operator == (const Point& b) const { return sgn(x - b.x) == 0 && sgn(y - b.y) == 0; }
    bool operator < (const Point& b) const {
        if (sgn(x - b.x) != 0) return x < b.x;
        return y < b.y;
    }
    Point operator - () const { return Point(-x, -y); }
    Point& operator += (const Point& b) { return *this = *this + b; }
    Point& operator -= (const Point& b) { return *this = *this - b; }
    Point& operator *= (ld b) { return *this = *this * b; }
    
    ld len2() const { return (*this) * (*this); }
    ld len() const { return sqrtl(len2()); }
    ld dis(const Point& b) const { return (*this - b).len(); }
    ld dis2(const Point& b) const { return (*this - b).len2(); }
    
    // 逆时针旋转 rad 弧度
    Point rot(ld rad) const {
        return Point(x * cosl(rad) - y * sinl(rad), x * sinl(rad) + y * cosl(rad));
    }
    // 逆时针旋转，已知 cos 和 sin
    Point rot(ld cosr, ld sinr) const {
        return Point(x * cosr - y * sinr, x * sinr + y * cosr);
    }
    Point rot90() const { return Point(-y, x); }
    
    // 向量夹角 [0, PI]
    ld ang(const Point& b) const {
        return acosl(max((ld)-1.0, min((ld)1.0, ((*this) * b) / (len() * b.len()))));
    }
    
    // to-left 测试: 1 左侧, 0 共线, -1 右侧
    int toleft(const Point& p) const {
        ld ans = (*this) % p;
        return (ans > EPS) - (ans < -EPS);
    }
    
    // 象限: 0原点 1x正 2第一象限 3y正 4第二象限 5x负 6第三象限 7y负 8第四象限
    int quad() const {
        if (sgn(x) == 0 && sgn(y) == 0) return 0;
        if (sgn(y) == 0) return x > 0 ? 1 : 5;
        if (sgn(x) == 0) return y > 0 ? 3 : 7;
        if (y > 0) return x > 0 ? 2 : 4;
        return x > 0 ? 8 : 6;
    }
};

// 两点距离
ld dis(const Point& a, const Point& b) {
    return sqrtl((a - b) * (a - b));
}

// 全平面极角排序
struct Argcmp {
    Point o;
    Argcmp(Point p = Point(0, 0)) : o(p) {}
    bool operator()(const Point& a, const Point& b) const {
        int qa = (a - o).quad(), qb = (b - o).quad();
        if (qa != qb) return qa < qb;
        ld t = (a - o) % (b - o);
        return t > EPS;
    }
};

直线与线段

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
struct Line {
    Point p, v;  // 点 + 方向向量
    Line() {}
    Line(Point a, Point b) : p(a), v(b) {}
    
    bool operator < (const Line& a) const {
        if (sgn(v % a.v) == 0 && v * a.v >= -EPS)
            return toleft(a.p) == -1;
        return Argcmp()(v, a.v);
    }
    
    bool is_connect(const Line& a) const {
        return sgn(v % a.v) != 0;
    }
    
    int toleft(const Point& a) const {
        return v.toleft(a - p);
    }
    
    ld dis(const Point& a) const {
        return fabs(v % (a - p)) / v.len();
    }
    
    Point proj(const Point& a) const {
        return p + v * ((v * (a - p)) / (v * v));
    }
    
    Point inter(const Line& a) const {
        assert(is_connect(a));
        return p + v * ((a.v % (p - a.p)) / (v % a.v));
    }
};

struct Segment {
    Point a, b;
    Segment() {}
    Segment(Point _a, Point _b) : a(_a), b(_b) {}
    
    // -1 端点, 0 不在, 1 严格在线段上
    int is_on(Point p) const {
        if (p == a || p == b) return -1;
        return (b - a).toleft(p - a) == 0 && (p - b) * (p - a) < -EPS;
    }
    
    // 线段与直线: -1 经过端点, 0 不相交, 1 严格相交
    int is_inter(const Line& l) const {
        if (l.toleft(a) == 0 || l.toleft(b) == 0) return -1;
        return l.toleft(a) != l.toleft(b);
    }
    
    // 两线段: -1 端点相交, 0 不相交, 1 严格相交
    int is_inter(const Segment& s) const {
        if (is_on(s.a) || is_on(s.b) || s.is_on(a) || s.is_on(b)) return -1;
        Line l(a, b - a), ls(s.a, s.b - s.a);
        return l.toleft(s.a) * l.toleft(s.b) == -1 && ls.toleft(a) * ls.toleft(b) == -1;
    }
    
    ld dis(const Point& p) const {
        if ((p - a) * (b - a) < -EPS || (p - b) * (a - b) < -EPS)
            return min(a.dis(p), b.dis(p));
        return Line(a, b - a).dis(p);
    }
    
    ld dis(const Segment& s) const {
        if (is_inter(s)) return 0;
        return min({dis(s.a), dis(s.b), s.dis(a), s.dis(b)});
    }
};

多边形与凸包

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
struct Polygon {
    vector<Point> point;
    Polygon() {}
    Polygon(vector<Point> p) : point(p) {}
    
    int nxt(int i) const { return i + 1 == (int)point.size() ? 0 : i + 1; }
    int pre(int i) const { return i == 0 ? (int)point.size() - 1 : i - 1; }
    
    // 回转数: {是否在边上, 回转数}
    pair<bool, int> winding(const Point& p) const {
        int cnt = 0;
        for (int i = 0; i < (int)point.size(); ++i) {
            Point u = point[i], v = point[nxt(i)];
            if (Segment(u, v).is_on(p)) return {true, 0};
            if (sgn(u.y - v.y) == 0) continue;
            Line uv(u, v - u);
            if (u.y < v.y - EPS && uv.toleft(p) <= 0) continue;
            if (u.y > v.y + EPS && uv.toleft(p) >= 0) continue;
            if (u.y < p.y - EPS && v.y >= p.y - EPS) cnt++;
            if (u.y >= p.y - EPS && v.y < p.y - EPS) cnt--;
        }
        return {false, cnt};
    }
    
    // 面积两倍，正表示逆时针
    ld area() const {
        ld sum = 0;
        for (int i = 0; i < (int)point.size(); ++i)
            sum += point[i] % point[nxt(i)];
        return sum;
    }
    
    ld circ() const {
        ld sum = 0;
        for (int i = 0; i < (int)point.size(); ++i)
            sum += point[i].dis(point[nxt(i)]);
        return sum;
    }
};

struct Convex : Polygon {
    Convex() {}
    Convex(vector<Point> p) : Polygon(p) {}
    
    // 闵可夫斯基和
    Convex operator + (const Convex& c) const {
        vector<Segment> e1(point.size()), e2(c.point.size()), edge(point.size() + c.point.size());
        vector<Point> res;
        res.reserve(point.size() + c.point.size());
        auto cmp = [](const Segment& u, const Segment& v) {
            return Argcmp()(u.b - u.a, v.b - v.a);
        };
        for (int i = 0; i < (int)point.size(); ++i)
            e1[i] = Segment(point[i], point[nxt(i)]);
        for (int i = 0; i < (int)c.point.size(); ++i)
            e2[i] = Segment(c.point[i], c.point[c.nxt(i)]);
        rotate(e1.begin(), min_element(e1.begin(), e1.end(), cmp), e1.end());
        rotate(e2.begin(), min_element(e2.begin(), e2.end(), cmp), e2.end());
        merge(e1.begin(), e1.end(), e2.begin(), e2.end(), edge.begin(), cmp);
        auto check = [](const vector<Point>& res, const Point& u) {
            Point back1 = res.back(), back2 = *prev(res.end(), 2);
            return (back1 - back2).toleft(u - back1) == 0 && (back1 - back2) * (u - back1) >= -EPS;
        };
        Point u = e1[0].a + e2[0].a;
        for (auto& v : edge) {
            while (res.size() > 1 && check(res, u)) res.pop_back();
            res.push_back(u);
            u = u + v.b - v.a;
        }
        if (res.size() > 1 && check(res, res[0])) res.pop_back();
        return Convex(res);
    }
    
    Convex operator - (const Convex& c) const {
        Convex nc = c;
        for (auto& p : nc.point) p = -p;
        return (*this) + nc;
    }
    
    // 旋转卡壳求凸包直径平方
    ld rotcaliper() const {
        if (point.size() <= 1) return 0;
        if (point.size() == 2) return point[0].dis2(point[1]);
        ld ret = 0;
        auto area = [](const Point& u, const Point& v, const Point& w) {
            return (w - u) % (w - v);
        };
        int n = point.size(), r = 1;
        for (int i = 0; i < n; ++i) {
            ret = max({ret, point[i].dis2(point[r]), point[nxt(i)].dis2(point[r])});
            while (area(point[i], point[nxt(i)], point[r]) <= area(point[i], point[nxt(i)], point[nxt(r)])) {
                r = nxt(r);
                ret = max({ret, point[i].dis2(point[r]), point[nxt(i)].dis2(point[r])});
            }
        }
        return ret;
    }
    
    // 判断点是否在凸多边形内: -1边上 0外 1内
    int is_in(const Point& p) const {
        int n = point.size();
        if (n == 1) return p == point[0] ? -1 : 0;
        if (n == 2) return Segment(point[0], point[1]).is_on(p) ? -1 : 0;
        if (p == point[0]) return -1;
        if ((point[1] - point[0]).toleft(p - point[0]) == -1) return 0;
        if ((point.back() - point[0]).toleft(p - point[0]) == 1) return 0;
        auto cmp = [&](const Point& u, const Point& v) {
            return (u - point[0]).toleft(v - point[0]) == 1;
        };
        int i = lower_bound(point.begin() + 1, point.end(), p, cmp) - point.begin();
        if (i == 1) return Segment(point[0], point[i]).is_on(p) ? -1 : 0;
        if (i == n - 1 && Segment(point[0], point[i]).is_on(p)) return -1;
        if (Segment(point[i - 1], point[i]).is_on(p)) return -1;
        return (point[i] - point[i - 1]).toleft(p - point[i - 1]) > 0;
    }
    
    // 凸多边形关于某一方向的极点，O(logn)
    template <typename F>
    int extreme(const F& dir) const {
        auto check = [&](int i) {
            return dir(point[i]).toleft(point[nxt(i)] - point[i]) >= 0;
        };
        auto dir0 = dir(point[0]);
        bool check0 = check(0);
        if (!check0 && check(point.size() - 1)) return 0;
        auto cmp = [&](const Point& v) {
            int vi = &v - point.data();
            if (vi == 0) return 1;
            bool checkv = check(vi);
            ld t = dir0 % (v - point[0]);
            if (vi == 1 && checkv == check0 && sgn(t) == 0) return 1;
            return checkv ^ (checkv == check0 && t <= 0);
        };
        return partition_point(point.begin(), point.end(), cmp) - point.begin();
    }
    
    // 凸多边形外一点的两条切线
    pair<int, int> tangent(const Point& a) const {
        if (point.size() == 1) return {0, 0};
        if (point.size() == 2) return {0, 1};
        int u = extreme([&](const Point& u) { return u - a; });
        int v = extreme([&](const Point& u) { return a - u; });
        return {u, v};
    }
    
    // 平行于给定直线的凸多边形切线
    pair<int, int> tangent(const Line& a) const {
        int u = extreme([&](const Point&) { return a.v; });
        int v = extreme([&](const Point&) { return -a.v; });
        return {u, v};
    }
};

// Andrew 算法求凸包，结果逆时针
Convex convexhull(vector<Point> p) {
    vector<Point> st;
    if (p.empty()) return Convex(st);
    if (p.size() == 1) return Convex(p);
    sort(p.begin(), p.end());
    auto check = [](const vector<Point>& st, const Point& u) {
        Point back1 = st.back(), back2 = *prev(st.end(), 2);
        return (back1 - back2).toleft(u - back1) <= 0;
    };
    for (const Point& u : p) {
        while (st.size() > 1 && check(st, u)) st.pop_back();
        st.push_back(u);
    }
    int k = st.size();
    p.pop_back();
    reverse(p.begin(), p.end());
    for (const Point& u : p) {
        while ((int)st.size() > k && check(st, u)) st.pop_back();
        st.push_back(u);
    }
    st.pop_back();
    return Convex(st);
}

// 半平面交，返回直线集合，lim 为边界范围
vector<Line> halfinter(vector<Line> l, ld lim = 1e18) {
    auto check = [](const Line& a, const Line& b, const Line& c) {
        return a.toleft(b.inter(c)) < 0;
    };
    l.push_back(Line(Point(-lim, 0), Point(0, -1)));
    l.push_back(Line(Point(0, -lim), Point(1, 0)));
    l.push_back(Line(Point(lim, 0), Point(0, 1)));
    l.push_back(Line(Point(0, lim), Point(-1, 0)));
    sort(l.begin(), l.end());
    deque<Line> q;
    for (int i = 0; i < (int)l.size(); ++i) {
        if (i > 0 && l[i - 1].v.toleft(l[i].v) == 0 && l[i - 1].v * l[i].v > EPS)
            continue;
        while (q.size() > 1 && check(l[i], q.back(), q[q.size() - 2]))
            q.pop_back();
        while (q.size() > 1 && check(l[i], q[0], q[1]))
            q.pop_front();
        if (!q.empty() && q.back().v.toleft(l[i].v) <= 0)
            return vector<Line>();
        q.push_back(l[i]);
    }
    while (q.size() > 1 && check(q[0], q.back(), q[q.size() - 2]))
        q.pop_back();
    while (q.size() > 1 && check(q.back(), q[0], q[1]))
        q.pop_front();
    return vector<Line>(q.begin(), q.end());
}

// 平面最近点对
ld closest_pair(vector<Point> points) {
    sort(points.begin(), points.end());
    auto cmpy = [](const Point& a, const Point& b) {
        if (sgn(a.y - b.y) == 0) return a.x < b.x - EPS;
        return a.y < b.y - EPS;
    };
    multiset<Point, decltype(cmpy)> s(cmpy);
    ld ans = INF;
    for (int i = 0, l = 0; i < (int)points.size(); ++i) {
        ld sqans = sqrtl(ans) + 1;
        while (l < i && points[i].x - points[l].x >= sqans)
            s.erase(s.find(points[l++]));
        for (auto it = s.lower_bound(Point(-INF, points[i].y - sqans));
             it != s.end() && it->y - points[i].y <= sqans; ++it)
            ans = min(ans, points[i].dis2(*it));
        s.insert(points[i]);
    }
    return ans;
}

圆相关

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
struct Circle {
    Point o;
    ld r;
    Circle() {}
    Circle(Point p, ld _r) : o(p), r(_r) {}
    
    bool operator == (const Circle& a) const {
        return o == a.o && sgn(r - a.r) == 0;
    }
    ld circ() const { return 2 * PI * r; }
    ld area() const { return PI * r * r; }
    
    // 点与圆: -1 圆上 0 圆外 1 圆内
    int is_in(const Point& p) const {
        ld d2 = p.dis2(o);
        if (sgn(d2 - r * r) == 0) return -1;
        return d2 < r * r - EPS;
    }
    
    // 直线与圆: 0 相离 1 相切 2 相交
    int relation(const Line& l) const {
        ld d = l.dis(o);
        if (d > r + EPS) return 0;
        if (sgn(d - r) == 0) return 1;
        return 2;
    }
    
    // 圆与圆: -1相同 0相离 1外切 2相交 3内切 4内含
    int relation(const Circle& a) const {
        if (*this == a) return -1;
        ld d = o.dis(a.o);
        if (d > r + a.r + EPS) return 0;
        if (sgn(d - r - a.r) == 0) return 1;
        if (sgn(d - fabs(r - a.r)) == 0) return 3;
        if (d < fabs(r - a.r) - EPS) return 4;
        return 2;
    }
    
    vector<Point> inter(const Line& l) const {
        ld d = l.dis(o);
        Point p = l.proj(o);
        int t = relation(l);
        if (t == 0) return vector<Point>();
        if (t == 1) return vector<Point>{p};
        ld k = sqrtl(r * r - d * d);
        Point dir = l.v / l.v.len();
        return vector<Point>{p - dir * k, p + dir * k};
    }
    
    vector<Point> inter(const Circle& a) const {
        ld d = o.dis(a.o);
        int t = relation(a);
        if (t == -1 || t == 0 || t == 4) return vector<Point>();
        Point e = (a.o - o) / (a.o - o).len() * r;
        if (t == 1 || t == 3) {
            if (r * r + d * d - a.r * a.r >= -EPS)
                return vector<Point>{o + e};
            return vector<Point>{o - e};
        }
        ld costh = (r * r + d * d - a.r * a.r) / (2 * r * d);
        ld sinth = sqrtl(1 - costh * costh);
        return vector<Point>{o + e.rot(costh, -sinth), o + e.rot(costh, sinth)};
    }
    
    ld interArea(const Circle& a) const {
        ld d = o.dis(a.o);
        int t = relation(a);
        if (t == -1) return area();
        if (t < 2) return 0;
        if (t > 2) return min(area(), a.area());
        ld costh1 = (r * r + d * d - a.r * a.r) / (2 * r * d);
        ld costh2 = (a.r * a.r + d * d - r * r) / (2 * a.r * d);
        ld sinth1 = sqrtl(1 - costh1 * costh1);
        ld sinth2 = sqrtl(1 - costh2 * costh2);
        ld th1 = acosl(costh1), th2 = acosl(costh2);
        return r * r * (th1 - costh1 * sinth1) + a.r * a.r * (th2 - costh2 * sinth2);
    }
    
    // 圆与多边形面积交
    ld interArea(const Polygon& plg) {
        auto cal = [&](const Point& a, const Point& b) {
            if (sgn((a - o) % (b - o)) == 0) return (ld)0;
            bool ina = is_in(a) == 1, inb = is_in(b) == 1;
            Line ab(a, b - a);
            if (ina && inb) return ((a - o) % (b - o)) / 2;
            if (ina && !inb) {
                auto pts = inter(ab);
                Point p = pts.size() == 1 ? pts[0] : pts[1];
                ld ans = ((a - o) % (p - o)) / 2;
                ld th = (p - o).ang(b - o);
                ld d = r * r * th / 2;
                return ((a - o) % (b - o)) > 0 ? ans + d : ans - d;
            }
            if (!ina && inb) {
                Point p = inter(ab)[0];
                ld ans = ((p - o) % (b - o)) / 2;
                ld th = (a - o).ang(p - o);
                ld d = r * r * th / 2;
                return ((a - o) % (b - o)) > 0 ? ans + d : ans - d;
            }
            auto pts = inter(ab);
            if (pts.size() == 2 && Segment(a, b).dis(o) <= r + EPS) {
                ld ans = ((pts[0] - o) % (pts[1] - o)) / 2;
                ld th1 = (a - o).ang(pts[0] - o);
                ld th2 = (b - o).ang(pts[1] - o);
                ld d1 = r * r * th1 / 2, d2 = r * r * th2 / 2;
                return ((a - o) % (b - o)) > 0 ? ans + d1 + d2 : ans - d1 - d2;
            }
            ld th = (a - o).ang(b - o);
            return ((a - o) % (b - o)) > 0 ? r * r * th / 2 : -r * r * th / 2;
        };
        ld ans = 0;
        for (int i = 0; i < (int)plg.point.size(); ++i) {
            Point a = plg.point[i], b = plg.point[plg.nxt(i)];
            ans += cal(a, b);
        }
        return ans;
    }
    
    // 圆外一点的切线
    vector<Line> tangent(const Point& p) const {
        int t = is_in(p);
        if (t == 1) return vector<Line>();
        if (t == -1) {
            Point v = Point(-(p - o).y, (p - o).x);
            return vector<Line>{Line(p, v)};
        }
        Point e = (p - o) / (p - o).len() * r;
        ld costh = r / o.dis(p);
        ld sinth = sqrtl(1 - costh * costh);
        Point t1 = o + e.rot(costh, -sinth);
        Point t2 = o + e.rot(costh, sinth);
        return vector<Line>{Line(p, t1 - p), Line(p, t2 - p)};
    }
    
    // 两圆公切线
    vector<Line> tangent(const Circle& a) const {
        int t = relation(a);
        vector<Line> lines;
        if (t == -1 || t == 4) return lines;
        if (t == 1 || t == 3) {
            Point p = inter(a)[0];
            Point v = Point(-(a.o - o).y, (a.o - o).x);
            lines.push_back(Line(p, v));
        }
        ld d = o.dis(a.o);
        Point e = (a.o - o) / (a.o - o).len();
        if (t <= 2) {
            ld costh = (r - a.r) / d;
            ld sinth = sqrtl(1 - costh * costh);
            Point d1 = e.rot(costh, -sinth), d2 = e.rot(costh, sinth);
            Point u1 = o + d1 * r, u2 = o + d2 * r;
            Point v1 = a.o + d1 * a.r, v2 = a.o + d2 * a.r;
            lines.push_back(Line(u1, v1 - u1));
            lines.push_back(Line(u2, v2 - u2));
        }
        if (t == 0) {
            ld costh = (r + a.r) / d;
            ld sinth = sqrtl(1 - costh * costh);
            Point d1 = e.rot(costh, -sinth), d2 = e.rot(costh, sinth);
            Point u1 = o + d1 * r, u2 = o + d2 * r;
            Point v1 = a.o - d1 * a.r, v2 = a.o - d2 * a.r;
            lines.push_back(Line(u1, v1 - u1));
            lines.push_back(Line(u2, v2 - u2));
        }
        return lines;
    }
};

常用几何算法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
// 三点外接圆
Circle CircumCircle(Point A, Point B, Point C) {
    Line L1((A + B) / 2, (B - A).rot90());
    Line L2((A + C) / 2, (C - A).rot90());
    if (!L1.is_connect(L2)) return Circle(Point(0, 0), INF);
    Point center = L1.inter(L2);
    return Circle(center, center.dis(A));
}

// 三点内切圆
Circle InCircle(Point A, Point B, Point C) {
    ld a = B.dis(C), b = A.dis(C), c = A.dis(B);
    ld area = fabs((B - A) % (C - A));
    Point center = (A * a + B * b + C * c) / (a + b + c);
    ld r = area / (a + b + c);
    return Circle(center, r);
}

// 最小圆覆盖
Circle min_cover_circle(vector<Point> ps) {
    int n = ps.size();
    random_shuffle(ps.begin(), ps.end());
    Circle cir(ps[0], 0);
    for (int i = 0; i < n; i++) {
        if (cir.is_in(ps[i]) != 1) continue;
        cir.o = ps[i], cir.r = 0;
        for (int j = 0; j < i; j++) {
            if (cir.is_in(ps[j]) != 1) continue;
            cir.o = (ps[i] + ps[j]) / 2;
            cir.r = ps[i].dis(ps[j]) / 2;
            for (int k = 0; k < j; k++) {
                if (cir.is_in(ps[k]) != 1) continue;
                cir = CircumCircle(ps[i], ps[j], ps[k]);
            }
        }
    }
    return cir;
}

// 阿波罗尼斯圆
Circle getApolloniusCircle(Point A, Point B, ld k) {
    assert(sgn(k - 1.0) != 0 && k > EPS);
    ld k2 = k * k;
    Point o = (B - A * k2) / (1.0 - k2);
    ld r = (k * A.dis(B)) / fabs(1.0 - k2);
    return Circle(o, r);
}

// 三角形面积两倍（带符号）
ld triangle_area2(const Point& a, const Point& b, const Point& c) {
    return (b - a) % (c - a);
}

bool in_triangle_closed(const Point& p, const Point& a, const Point& b, const Point& c) {
    return triangle_area2(a, b, p) >= -EPS &&
           triangle_area2(b, c, p) >= -EPS &&
           triangle_area2(c, a, p) >= -EPS;
}

// 多边形三角剖分（待实现）
// vector<Triangle> polygon_triangulation(const Polygon& plg) {
//     // Triangle 结构需要额外定义
// }

参考资料

  • 本文档内容集成自 常用板子 模板
  • 更多算法模板请参考主模板文档