Visits: 564

一直莫名不想碰的一个专题

基础知识

点、向量模板

const ld eps = 1e-10L, PI = acos(-1);

struct P {
    ld x, y;
    inline P() {}
    inline P(ld x, ld y) : x(x), y(y) {}
    inline P &operator += (P o) { return x += o.x, y += o.y, *this; }
    inline P &operator -= (P o) { return x -= o.x, y -= o.y, *this; }
    inline P &operator *= (ld o) { return x *= o, y *= o, *this; }
    inline P &operator /= (ld o) { return x /= o, y /= o, *this; }
    inline friend P operator + (P a, P b) { return a += b; }
    inline friend P operator - (P a, P b) { return a -= b; }
    inline friend P operator * (P a, ld b) { return a *= b; }
    inline friend P operator / (P a, ld b) { return a /= b; }
    inline friend bool operator < (P a, P b) { return fabs(a.x - b.x) < eps ? a.y < b.y : a.x < b.x; }
    inline friend bool operator == (P a, P b) { return !(a < b) && !(b < a); }
    inline friend ld operator * (P a, P b) { return a.x * b.x + a.y * b.y; }
    inline friend ld operator % (P a, P b) { return a.x * b.y - a.y * b.x; }
    inline friend ld operator ^ (P a, P b) { return a * b / a.l() / b.l(); }
    inline ld a() { return atan2(y, x); }
    inline ld l() { return sqrt(*this * *this); }
    inline void r(ld o) { ld s = sin(o), c = cos(o), X = x * c - y * s, Y = x * s + y * c; x = X, y = Y; }
};

向量点积

上面模板中的 * 运算。

两个向量 $p_1,p_2$ 的点积定义为 $p_1,p_2$ 的模长乘它们夹角的 $\cos$，可以解释为 $p_1,p_2$ 做的功的大小。

$\cos$ 是偶函数，因此点积满足交换律。

$$p_1 * p_2 = x_1x_2 + y_1y_2 = p_2 * p_1$$

向量叉积

上面模板中的 % 运算。

两个向量 $p_1,p_2$ 的叉积可以解释为 $p_1,p_2$ 叉出来的平行四边形的有向面积。

叉积没有交换律，但交换后为原来的相反数。

$$p_1 \times p_2 = x_1y_2 – x_2y_1 = -p_2 \times p_1$$

向量夹角

上面模板中的 ^ 运算。

向量极角

上面模板中的 a 运算。

向量模长

上面模板中的 l 运算。

向量旋转

上面模板中的 r 运算。

向量绕起点逆时针旋转角度 $a$，公式为 $X = x \cos(a) – y \sin(a), Y = x \sin(a) + y \cos(a)$，本质上是个矩阵乘法。

$$
\begin{bmatrix} X & Y
\end{bmatrix} =
\begin{bmatrix} x & y
\end{bmatrix} \times
\begin{bmatrix}\text{cos}(a) & \text{sin}(a)\\ -\text{sin}(a) & \text{cos}(a)\end{bmatrix}
$$

直线模板

struct L {
    P a, b;
    inline L() {}
    inline L(P a, P b) : a(a), b(b) {}
    inline friend P operator * (L a, L b) { return a.a + a.b * (b.b % (a.a - b.a) / (a.b % b.b)); }
    inline friend ld operator / (L a, P b) { return fabs(a.b % (b - a.a)) / a.b.l(); }
    inline friend bool operator < (L a, L b) { return fabs(a.b.a() - b.b.a()) < eps ? (b.a - a.a) % a.b < 0 : a.b.a() < b.b.a(); }
};

如果被卡常，可以把 ld 换成 double，同时提前存下 b.a() 的值。

求两条直线的交点

上面模板中的 * 运算。

设直线 $a$ 的关键点 $a.a = P$，方向向量 $a.b = V$，直线 $b$ 的关键点 $b.a = Q$，方向向量 $b.b = W$。

所以交点 $(x,y)$ 满足 $x = P_x + t_1V_x = Q_x + t_2W_x, y = P_y + t_1V_y = Q_y + t_2W_y$。

解得 $t_1 = \frac{(P_y – Q_y)W_x – (P_x – Q_x)W_y}{V_xW_y – V_yW_x} = (W \times (P – Q) / (V \times W)) = (b.b \times (a.a – b.a) / (a.b \times b.b))$。

然后…就做完了…

顺便一提，线段交点就把直线交点求出来然后判一下就完了…

点到直线的距离

上面模板中的 / 运算。

用叉积（平行四边形面积）除以低即可。

欧拉定理

平面图满足：点数 $+$ 面数 $=$ 边数 $+ 2$。

判断一个点是否在任意多边形内部

从这个点出发引一条射线，如果这条射线与多边形有奇数个交点则在内部，否则在外部。

为了避免射线与多边形上的点重合，可以将射线的斜率设为无理数（随便乱打一个数）。

三角剖分求多边形面积

把多边形三角剖分，分别计算有向面积，然后求和。

二维凸包

二维上把若干给定点包围在内部的、面积最小的凸多边形。

先排序，然后单调栈维护一下就好了，时间复杂度 $\mathcal O(n \log n)$。

【模板】P2742 【模板】二维凸包 / [USACO5.1]圈奶牛Fencing the Cows

inline vector<P> convex_hull(vector<P> a) {
    sort(a.begin(), a.end());
    int n = unique(a.begin(), a.end()) - a.begin();
    vector<P> p;
    for (int i = 0; i < n; i++) {
        while (p.size() > 1u &&
            (p.back() - p[p.size()-2]) % (a[i] - p[p.size()-2]) < eps)
                p.pop_back();
        p.pb(a[i]);
    }
    ui m = p.size();
    for (int i = n - 2; ~i; i--) {
        while (p.size() > m &&
            (p.back() - p[p.size()-2]) % (a[i] - p[p.size()-2]) < eps)
                p.pop_back();
        p.pb(a[i]);
    }
    p.pop_back();
    return p;
}

int main() {
    int n;
    scanf("%d", &n);
    vector<P> a;
    ld x, y;
    for (int i = 1; i <= n; i++)
        scanf("%Lf %Lf", &x, &y), a.pb(P(x, y));
    a = convex_hull(a);
    a.pb(a[0]);
    ld ans = 0;
    for (ui i = 1; i < a.size(); i++) ans += (a[i] - a[i-1]).l();
    printf("%.2Lf\n", ans);
    return 0;
}

旋转卡壳

用于求平面最远点对。

首先求出凸包，目的是把所有点按逆时针排序。然后枚举逆时针边，定义 $T$ 为当前的对踵点，若 $T$ 逆时针变动能增大到当前边的距离则改变 $T$，直到继续改变会导致距离减小为止。每次求出边（设为 $AB$）的对踵点 $T_i$ 后，分别用 $|AT_i|,|BT_i|$ 来更新答案，即凸多边形的直径。

（对踵点：如果过凸多边形上两点作一对平行线，使得整个多边形都在这两条线之间，那么这两个点被称为一对对踵点。）

时间复杂度 $\mathcal O(n \log n)$。

【模板】P1452 Beauty Contest

const int N = 5e4 + 7;
int n, m;
P a[N], p[N];

inline ld S(P x, P a, P b) { return Abs((a - x) % (b - x) / (a - b).l()); }

int main() {
    ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
    cin >> n;
    for (int i = 1; i <= n; i++) cin >> a[i].x >> a[i].y;
    sort(a + 1, a + n + 1);
    n = unique(a + 1, a + n + 1) - (a + 1);
    for (int i = 1; i <= n; i++) {
        while (m > 1 && (p[m] - p[m-1]) % (a[i] - p[m-1]) <= 0) --m;
        p[++m] = a[i];
    }
    int o = m;
    for (int i = n - 1; i; i--) {
        while (m > o && (p[m] - p[m-1]) % (a[i] - p[m-1]) <= 0) --m;
        p[++m] = a[i];
    }
    if (m < 3) return cout << 0 << endl, 0;
    if (m == 3) {
        ld ans = (p[1] - p[2]).l();
        cout << (ll)(ans * ans + 0.5) << endl;
        return 0;
    }
    int t = 1;
    ld ans = 0;
    for (int i = 1; i < m; i++) {
        while (S(p[t], p[i], p[i+1]) <= S(p[t+1], p[i], p[i+1])) t = t % (m - 1) + 1;
        ans = Max(ans, Max((p[t] - p[i]).l(), (p[t] - p[i+1]).l()));
    }
    cout << (ll)(ans * ans + 0.5) << endl;
    return 0;
}

注意当凸包退化成一个点或一条线时需要特判。

半平面交

给出若干个半平面（有向直线的左侧），求它们的交。

做法与二维凸包比较像。先极角排序，然后单调队列维护即可，时间复杂度 $\mathcal O(n \log n)$。

【模板】P4196 [CQOI2006]凸多边形

半平面交存在三种情况：无解，有解且有界，有解但无界。

下面的 half_plane 中，如果无解，则返回的 deque<P> 大小为 $0$，否则为半平面交。

本题中不可能出现有解但无界的情况，如果可能出现，我们需要在最外面套上一个巨大的正方形转化成有界的情况，最后如果有解，但有边在外层的正方形上，则说明无界。

inline deque<P> half_plane(vector<L> l) {
    sort(l.begin(), l.end());
    deque<P> r;
    deque<L> q;
    for (L o : l) {
        while (r.size() && (r.back() - o.a) % o.b > 0)
            q.pop_back(), r.pop_back();
        while (r.size() && (r.front() - o.a) % o.b > 0)
            q.pop_front(), r.pop_front();
        if (q.size() && fabs(o.b % q.back().b) < eps) {
            if (fabs(o.b.a() - q.back().b.a()) > eps)
                return r.clear(), r;
            q.pop_back();
            if (r.size()) r.pop_back();
        }
        if (q.size()) r.pb(q.back() * o);
        q.pb(o);
    }
    while (r.size() && (r.back() - q.front().a) % q.front().b > 0)
        q.pop_back(), r.pop_back();
    if (q.size() < 3u) r.clear();
    else r.pb(q.front() * q.back());
    return r;
}

int main() {
    int n;
    vector<P> p;
    vector<L> l;
    rd(n);
    for (int i = 1, m, x, y; i <= n; i++) {
        rd(m), p.clear();
        while (m--) rd(x), rd(y), p.pb(P(x, y));
        p.pb(p[0]);
        for (ui i = 1; i < p.size(); i++) l.pb(L(p[i-1], p[i] - p[i-1]));
    }
    deque<P> r = half_plane(l);
    if (!r.size()) return prints("0.000"), 0;
    ld ans = r.back() % r.front();
    while (r.size() > 1u) {
        P o = r.back();
        r.pop_back();
        ans += r.back() % o;
    }
    cout << fixed << setprecision(3) << ans / 2 << endl;
    return 0;
}

扫描线

这一一个跟计算几何没关系却又有关系的算法…

线段树维护，由于扫描线的特殊性质可以不用下传标记。

【模板】P5490 【模板】扫描线

#define ls (p << 1)
#define rs (ls | 1)
#define mid ((t[p].l + t[p].r) >> 1)
#define ll long long

const int N = 1e5 + 6;
int n;
struct P {
    int x, y, z, k;
    inline bool operator < (const P &o) const {
        return x < o.x;
    }
} a[N<<1];
int raw[N<<1];
map<int, int> val;
struct T {
    int l, r, cnt, len;
} t[N<<3];
ll ans;

void build(int p, int l, int r) {
    t[p].l = l, t[p].r = r;
    if (l == r) return;
    build(ls, l, mid), build(rs, mid + 1, r);
}

void change(int p, int l, int r, int k) {
    if (l <= t[p].l && r >= t[p].r) return t[p].len = ((t[p].cnt += k) ? raw[t[p].r+1] - raw[t[p].l] : (t[p].l == t[p].r ? 0 : t[ls].len + t[rs].len)), void();
    if (l <= mid) change(ls, l, r, k);
    if (r > mid) change(rs, l, r, k);
    t[p].len = (t[p].cnt ? raw[t[p].r+1] - raw[t[p].l] : t[ls].len + t[rs].len);
}

int main() {
    rd(n);
    for (int i = 1; i <= n; i++) {
        int k = i << 1, y, z;
        rd(a[k-1].x), rd(y), rd(a[k].x), rd(z);
        raw[k-1] = a[k-1].y = a[k].y = y;
        raw[k] = a[k-1].z = a[k].z = z;
        a[k-1].k = 1, a[k].k = -1;
    }
    n <<= 1;
    sort(raw + 1, raw + n + 1);
    int m = unique(raw + 1, raw + n + 1) - (raw + 1);
    for (int i = 1; i <= m; i++) val[raw[i]] = i;
    sort(a + 1, a + n + 1);
    build(1, 1, m - 1);
    for (int i = 1; i < n; i++) {
        int y = val[a[i].y], z = val[a[i].z] - 1;
        change(1, y, z, a[i].k);
        ans += 1ll * t[1].len * (a[i+1].x - a[i].x);
    }
    print(ans);
    return 0;
}

曼哈顿距离和切比雪夫距离

设 $A(x_1, y_1), B(x_2, y_2)$。

• 曼哈顿距离：$|x_{1}-x_{2}|+|y_{1}-y_{2}|$。
• 切比雪夫距离：$\max(|x_{1}-x_{2}|,|y_{1}-y_{2}|)$。

相互转化

• 曼哈顿距离 $\to$ 切比雪夫距离：$(x, y) \to (x + y, x – y)$。
• 切比雪夫距离 $\to$ 曼哈顿距离：$(x, y) \to (\frac{x+y}{2},\frac{x-y}{2})$。

【例题】P3964 [TJOI2013]松鼠聚会

const int N = 1e5 + 7;
int n, x[N], y[N], p[N], q[N];
ll sx[N], sy[N], ans = 4e18;

int main() {
    rd(n);
    for (int i = 1, a, b; i <= n; i++) rd(a), rd(b), x[i] = p[i] = a + b, y[i] = q[i] = a - b;
    sort(p + 1, p + n + 1), sort(q + 1, q + n + 1);
    for (int i = 1; i <= n; i++) sx[i] = sx[i-1] + p[i], sy[i] = sy[i-1] + q[i];
    for (int i = 1; i <= n; i++) {
        int px = lower_bound(p + 1, p + n + 1, x[i]) - p, py = lower_bound(q + 1, q + n + 1, y[i]) - q;
        ans = min(ans, 1ll * px * x[i] - sx[px] + sx[n] - sx[px] - 1ll * (n - px) * x[i] + 1ll * py * y[i] - sy[py] + sy[n] - sy[py] - 1ll * (n - py) * y[i]);
    }
    print(ans >> 1);
    return 0;
}

参考资料

• 刘汝佳 & 陈锋 《算法竞赛进阶指南·入门经典》
• 李煜东 《算法竞赛进阶指南》
• wjyyy 《计算几何 学习笔记合集【凸包】【旋转卡壳】【半平面交】》
• 自为风月马前卒 曼哈顿距离与切比雪夫距离及其相互转化
• mcfx 《计算几何》课件
• 吴清月 《计算几何》课件