「ACM-ICPC 2018 南京站网络赛 D 题」Jerome's House - 半平面交 + 三分

建议访问原文出处，获得更佳浏览体验。
原文出处：https://hyp1231.github.io/2018/09/04/20180904-2018nanjing-online-d/

题意

给一个有 $n$ 个点的凸多边形，让你把三个半径为 $r$ 的圆放入多边形内，要求圆不能覆盖住多边形的边（但圆与圆可以相互覆盖）。求 $|x_1y_2-x_1y_3 + x_2y_3 - x_2y_1+x_3y_1 - x_3y_2|$ 的最大值（其中 $(x_i, y_i)$ 为圆 $i$ 的圆心坐标）。

链接

D. Jerome’s House

题解

首先我们分析题目要求的目标函数，实际就是：

$$\left| \begin{array}{ccc} x_1 & y_1 & 1 \\ x_2 & y_2 & 1 \\ x_3 & y_3 & 1 \end{array} \right| = 2 \times S_{\Delta}$$
也就是说题目希望使三个圆心构成的三角形的面积最大。

现在我们考察三个圆心的活动范围。由于圆不能覆盖住多边形的边，所以最大值的情况只能是和边相切，也就是圆心活动的范围，距离边的距离不小于 $r$。因此我们考虑直接构造一个小多边形，使得新多边形的边和对应大多边形的边距离为 $r$。我们直接将原多边形的各个边所在的有向直线向多边形内侧移动 $r$，然后再对移动后的有向直线做半平面交，求得小多边形。

现在考虑何时三角形面积最大。显然当三个点都为小多边形顶点时，三角形面积最大。但直接暴力枚举的时间复杂度为 $O(n^3)$，对于 $n = 1000$ 我们无法接受。因此我们只枚举两个顶点，而这两个顶点把小多边形的顶点分为两部分。

我们只考虑其中一部分，从一端点到另一端点，三角形面积先变大后变小，凸性函数我们可以三分求得最值。但假如我们目前枚举的两个顶点是 $l$ 和 $r$，那我们知道可以在 $[l,r]$ 上做三分，但其余部分 $[0, l]$ 和 $[r, n -1 ]$ 怎么办呢？

这里用到一个小技巧，把原顶点序列 $[p_1, p_2, \ldots, p_m]$ 复制一份变成 $[p_1, p_2, \ldots, p_m, p_{m+1}, \ldots, p_{2m}]$（其中 $p_i = p_{i + m}$）。这样剩下的部分，$[r, l + m]$ 上的三角形面积也是个凸函数，这样就又可以三分啦！

时间复杂度，半平面交 $O(n\log n)$；查找最大面积时，暴力枚举 $O(n^2)$，三分 $O(\log_3n)$。总时间复杂度 $O(n\log n + n^2 \log_3n)$。

代码

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <vector>
#include <cmath>
using namespace std;

const double EPS = 1e-6;         //精度系数

struct Point {
    double x, y;
    Point(double x = 0, double y = 0) :x(x), y(y) {}
    const bool operator < (Point A)const {
        return x == A.x ? y < A.y : x < A.x;
    }
};  //点的定义

typedef Point Vector;   //向量的定义

Vector operator + (Vector A, Vector B) { return Vector(A.x + B.x, A.y + B.y); }    //向量加法
Vector operator - (Vector A, Vector B) { return Vector(A.x - B.x, A.y - B.y); }    //向量减法
Vector operator * (Vector A, double p) { return Vector(A.x*p, A.y*p); }        //向量数乘

inline int dcmp(double x) {
    if (fabs(x) < EPS)return 0; else return x < 0 ? -1 : 1;
}   //与0的关系

inline double Dot(Vector A, Vector B) { return A.x * B.x + A.y * B.y; }    //向量点乘
inline double Length(Vector A) { return sqrt(Dot(A, A)); }     //向量长度
inline double Cross(Vector A, Vector B) { return A.x*B.y - A.y*B.x; }    //向量叉乘

struct Line {
    Point p;
    Vector v;
    double ang;
    Line(){}
    Line(Point p, Vector v) :p(p), v(v) { ang = atan2(v.y, v.x); }
    Point point(double t) { return p + v*t; }    //给t值返回点坐标
    bool operator < (const Line& L) const {
        return ang < L.ang;
    }
};  //直线

// 点 p 在有向直线 L 的左边（线上不算）
inline bool OnLeft(Line L, Point p) {
    return dcmp(Cross(L.v, p - L.p)) > 0;
}

// 二直线交点。假设交点唯一存在
inline Point GetIntersection(Line a, Line b) {
    Vector u = a.p - b.p;
    double t = Cross(b.v, u) / Cross(a.v, b.v);
    return a.p + a.v * t;
}

// 半平面交主过程
inline int HalfplaneIntersection(Line *L, int n, Point *poly) {
    sort(L, L + n);// 按极角排序

    int first, last;                // 双端队列的第一个元素和最后一个元素的下标
    Point *p = new Point[n];        // p[i] 为 q[i] 和 q[i + 1] 的交点
    Line *q = new Line[n];          // 双端队列
    q[first = last = 0] = L[0];     // 双端队列初始化为只有一个半平面 L[0]
    for (int i = 1; i < n; ++i) {
        while (first < last && !OnLeft(L[i], p[last - 1])) --last;
        while (first < last && !OnLeft(L[i], p[first])) ++first;
        q[++last] = L[i];
        if (fabs(Cross(q[last].v, q[last - 1].v)) < EPS) {
        // 两向量平行且同向
            --last;
            if (OnLeft(q[last], L[i].p)) q[last] = L[i];
        }
        if (first < last) p[last - 1] = GetIntersection(q[last - 1], q[last]);
    }
    while (first < last && !OnLeft(q[first], p[last - 1])) --last;
    // 删除无用平面
    if (last - first <= 1) return 0;                // 空集
    p[last] = GetIntersection(q[last], q[first]);   // 计算首尾两个半平面的交点

    // 从双端队列复制到输出中
    int m = 0;
    for (int i = first; i <= last; ++i) poly[m++] = p[i];

    delete [] p; delete [] q;
    return m;
}

inline Vector Normal(Vector A) {
    double L = Length(A);
    return Vector(-A.y / L, A.x / L);
}   //单位法线，左转90°化为单一长度，确保A不是零向量

const int N = 1024;

Point input[N], poly[N << 1];
Line lines[N];

int L, R;

inline double f(int x) {
    Point &A = poly[L], &B = poly[R], &C = poly[x];
    return fabs(A.x * B.y - A.x * C.y + B.x * C.y - B.x * A.y + C.x * A.y - C.x * B.y);
}

inline double tri_divide(int l, int r) {
    while (l < r - 1) {
        int mid = (l + r) >> 1;
        int mmid = (mid + r) >> 1;
        if (dcmp(f(mid) - f(mmid)) > 0) {
            r = mmid;
        } else {
            l = mid;
        }
    }
    return std::max(f(l), f(r));
}

int main() {
    int T;
    scanf("%d", &T);
    while (T--) {
        int n;
        double ra;
        scanf("%d%lf", &n, &ra);
        for (int i = 0; i < n; ++i) {
            scanf("%lf%lf", &input[i].x, &input[i].y);
        }
        for (int i = 0; i < n; ++i) {
            lines[i] = Line(input[(i + 1) % n], input[i] - input[(i + 1) % n]);
        }   // 逆时针构造半平面在左侧的有向直线
        for (int i = 0; i < n; ++i) {
            lines[i].p = lines[i].p + Normal(lines[i].v) * ra;
        }   // 半平面向多边形内缩 r
        int m = HalfplaneIntersection(lines, n, poly);
            // 半平面交求缩过后的多边形
        for (int i = m; i < 2 * m; ++i) {
            poly[i] = poly[i - m];
        }   // 小技巧，复制一遍方便三分
        // 暴力枚举凸包上的两个点，三分求第三个点使面积最大
        double best = 0;
        for (int l = 0; l < m - 1; ++l) {
            for (int r = l + 1; r < m; ++r) {
                // 第一次三分，[l, r]
                L = l, R = r;
                double tmp = tri_divide(l, r);
                // 第二次三分，[0, l] + [r, n - 1]
                // 由于我们复制了一份多边形顶点序列，上述范围等价于 [r, l + m]
                L = r, R = l + m;
                tmp = std::max(tmp, tri_divide(l, r));
                best = std::max(best, tmp);
            }
        }
        printf("%.6f\n", best);
    }

    return 0;
}