关于812题的 \(O(n)\) 算法的简单思考
因为今天的题目很有意思所以特别想跟大家分享一下。
812. 最大三角形面积
一开始我想到了凸包，然后想到凸包后可以采用 \(O(n^2)\) 的渐进算法算出最大面积。但是灵神的回答中提到了一篇论文！
Maximal Area Triangles in a Convex Polygon
作者声称我们可以通过 \(O(n)\) 的算法复杂度找到这个最大值。
于是我就写出了下面的这一大段💩山代码
class Solution {
public:double largestTriangleArea(vector<vector<int>> &points) {int n = points.size();int *stack = new int[n + 2];memset(stack, -1, sizeof(int) * (n + 2));bool *used = new bool[n];memset(used, false, sizeof(bool) * n);ranges::sort(points, {}, [](const vector<int> &p) -> tuple<const int &, const int &> {return tie(p[0], p[1]);});// lower convex hull.int stack_pt = 0;for (int index = 0; index < (int)points.size(); index++) {int x = points[index][0], y = points[index][1];// When \vec{StSt-1 x StP <= 0} and more than 2 elements in stack: popwhile (stack_pt + 1 > 2) {int top1 = stack[stack_pt], top2 = stack[stack_pt - 1];int x1 = points[top1][0], y1 = points[top1][1];int x2 = points[top2][0], y2 = points[top2][1];if ((x1 - x2) * (y - y1) - (y1 - y2) * (x - x1) >= 0) {used[stack[stack_pt]] = false;stack_pt--;continue;} else {break;}}if (used[index] == false) {stack[++stack_pt] = index;used[index] = true;continue;}}// upper convex hull.int lower_size = stack_pt;for (int index = (int)points.size() - 1; index >= 0; index--) {int x = points[index][0], y = points[index][1];// When \vec{StSt-1 x StP <= 0} and more than 2 elements in stack: popwhile (stack_pt > lower_size) {int top1 = stack[stack_pt], top2 = stack[stack_pt - 1];int x1 = points[top1][0], y1 = points[top1][1];int x2 = points[top2][0], y2 = points[top2][1];if ((x1 - x2) * (y - y1) - (y1 - y2) * (x - x1) >= 0) {used[stack[stack_pt]] = false;stack_pt--;continue;} else {break;}}if (used[index] == false) {stack[++stack_pt] = index;used[index] = true;continue;}}// Calculate the triangle by calculate the convex hull point.auto S = [](int x1, int y1, int x2, int y2, int x3, int y3) -> double {return 0.5 * abs(x1 * y2 + x2 * y3 + x3 * y1 - x1 * y3 - x2 * y1 - x3 * y2);};auto dist = [&](int a, int b, int c) {int xa = points[stack[a]][0], ya = points[stack[a]][1];int xb = points[stack[b]][0], yb = points[stack[b]][1];int xc = points[stack[c]][0], yc = points[stack[c]][1];return (yb - ya) * xc + (xa - xb) * yc;};auto area = [&](int a, int b, int c) -> double {int x1 = points[stack[a]][0], y1 = points[stack[a]][1];int x2 = points[stack[b]][0], y2 = points[stack[b]][1];int x3 = points[stack[c]][0], y3 = points[stack[c]][1];return 0.5 * abs(x1 * y2 + x2 * y3 + x3 * y1 - x1 * y3 - x2 * y1 - x3 * y2);};double ans = 0;stack[0] = stack[stack_pt], stack[stack_pt + 1] = stack[1];if (stack_pt == 3) {int x1 = points[stack[0]][0], y1 = points[stack[0]][1];int x2 = points[stack[1]][0], y2 = points[stack[1]][1];int x3 = points[stack[2]][0], y3 = points[stack[2]][1];return S(x1, y1, x2, y2, x3, y3);}// Find one 3 stable.auto findOne3Stable = [&]() -> tuple<int, int, int> {double max_S = 0.0;int a = 1, b = 2, c = 3;int C = c;for (int B = b; B < stack_pt; B++) {if (C == B) C++;while (C < stack_pt && dist(1, B, C + 1) > dist(1, B, C) + 1e-6) {C++;}double new_area = area(1, B, C);if (new_area > max_S) {max_S = new_area;a = 1;b = B;c = C;}}// Step 2: find one 2-stable.while (dist(b, c, a + 1) > dist(b, c, a) + 1e-6) {a = a % stack_pt + 1;bool flag = true;while (flag == true) {flag = false;while (dist(c, a, b + 1) > dist(c, a, b) + 1e-6) {b = b % stack_pt + 1;flag = true;}while (dist(a, b, c + 1) > dist(a, b, c) + 1e-6) {c = c % stack_pt + 1;flag = true;}}}// Step 3: find one 3-stable.while (dist(b, c, a - 1) > dist(b, c, a) + 1e-6) {a--;if (a == 0) a = stack_pt;bool flag = true;while (flag == true) {flag = false;while (dist(c, a, b - 1) > dist(c, a, b) + 1e-6) {b--;if (b == 0) b = stack_pt;flag = true;}while (dist(a, b, c - 1) > dist(a, b, c) + 1e-6) {c--;if (c == 0) c = stack_pt;flag = true;}}}return {a, b, c};};auto RotateAndKill = [&](int a, int b, int c) -> double {int ans_a = 1, ans_b = 2, ans_c = 3;double max_a = 0.0;int r = a, t = c;// Section 4:while (b != t || c != r) {// repeat 1. find A.while (dist(b, c, a + 1) > dist(b, c, a)) {a = a % stack_pt + 1;}double new_area = area(a, b, c);// output the 3-stable.if (new_area > max_a) {max_a = new_area;ans_a = a, ans_b = b, ans_c = c;}// Repeat 2: Calculate Ib,c.double dx1 = points[stack[b]][0] - points[stack[b + 1]][0], dy1 = points[stack[b + 1]][1] - points[stack[b]][1];double dx2 = points[stack[c]][0] - points[stack[c + 1]][0], dy2 = points[stack[c + 1]][1] - points[stack[c]][1];double delta = dy1 * dx2 - dy2 * dx1;if (delta >= -1e-9)b = b % stack_pt + 1;else {double U = dy1 * points[stack[c]][0] + dx1 * points[stack[c]][1];double V = dy2 * points[stack[b]][0] + dx2 * points[stack[b]][1];double Ix = (U * dx2 - V * dx1) / delta;double Iy = (V * dy1 - U * dy2) / delta;// (yc - yb) * xa + (xb - xc) * ya;if (dist(b, c, a) > (points[stack[c]][1] - points[stack[b]][1]) * Ix + (points[stack[b]][0] - points[stack[c]][0]) * Iy) {c = c % stack_pt + 1;} else {b = b % stack_pt + 1;}}}a = ans_a, b = ans_b, c = ans_c;return max_a;};// Step 1:auto [a, b, c] = findOne3Stable();ans = RotateAndKill(a, b, c);return ans;}
};
我们姑且不管数学方面是否正确，我们来尝试部署一下这个算法。
前置知识:
- stable：对于一个已经提供了凸多边形点集的三角形子集 ABC，A 是 stable 的，当且仅当在固定了 BC 后，A 距离 BC 的距离是最大的。
- 我们所有的点将是顺时针的。即凸包上的点按照顺时针进行编号。
接下来我们来看算法是如何实现的。
- 找到一个 3-stable 三角形。也就是说我们首先要找到一个三角形 ABC，满足 A，B，C 都是 stable 的条件。
- 固定一个 A，然后寻找 B，C。这样我们就可以得到一个以 A 为根，B，C stable 的三角形。我们通过枚举 B，判断 C 是否是 stable 的。因为每次枚举 B，我们无需将 C 从头开始枚举，只需要从上一次的点继续向下枚举。因此摊还分析为 O(1) 平均操作（这一部分可以在官解上看到相同的做法）。对于 B 的枚举就是 \(O(n)\)。于是我们有：
double max_S = 0.0;int a = 1, b = 2, c = 3;int C = c;for (int B = b; B < stack_pt; B++) {if (C == B) C++;while (C < stack_pt && dist(1, B, C + 1) > dist(1, B, C) + 1e-6) {C++;}double new_area = area(1, B, C);if (new_area > max_S) {max_S = new_area;a = 1;b = B;c = C;}}
- 接着，作者假设了这个 A 并不是一个 stable 点。因此需要继续进行枚举，不断枚举根节点 A 来判断这个点是不是 stable 的。这里用了一个小 trick。我们的 stack 记录在 [1...stack_pt] 中，向外延拓两点：
stack[0] = stack[stack_pt]以及stack[stack_pt+1]=stack[1]来规避掉错误的地址访问。这样我们可以尝试讨论 A 使它成为 stable 点。我们从两个方向来分别判断它是不是稳定的- 顺时针：要求 A 顺时针方向的下一个节点 A+1到达 BC 的距离最大。
- 逆时针：要求A 逆时针方向的下一个节点 A-1到达 BC 的距离最大。
这里距离可以直接采用化简后的叉乘公式。我们可以回顾一下。这里我们统一为顺时针。固定住 AB，我们有：
这样可以让数值更加稳定。
// Step 2: find one 2-stable.while (dist(b, c, a + 1) > dist(b, c, a) + 1e-6) {a = a % stack_pt + 1;bool flag = true;while (flag == true) {flag = false;while (dist(c, a, b + 1) > dist(c, a, b) + 1e-6) {b = b % stack_pt + 1;flag = true;}while (dist(a, b, c + 1) > dist(a, b, c) + 1e-6) {c = c % stack_pt + 1;flag = true;}}}// Step 3: find one 3-stable.while (dist(b, c, a - 1) > dist(b, c, a) + 1e-6) {a--;if (a == 0) a = stack_pt;bool flag = true;while (flag == true) {flag = false;while (dist(c, a, b - 1) > dist(c, a, b) + 1e-6) {b--;if (b == 0) b = stack_pt;flag = true;}while (dist(a, b, c - 1) > dist(a, b, c) + 1e-6) {c--;if (c == 0) c = stack_pt;flag = true;}}}
这样就找到一个 3-Stable 了。
接下来进入到寻找所有的 stable 三角形。作者的做法是：
最关键的部分是\(I_{b,c}\)。我们发现 \(I_{b,c}\) 的选取方式是如下图所示的：
- 对于 C 点，做一条射线与 B 点和 B 的顺时针下一个点 B+1 的方向相反。
- 对于 B 点，做一条射线与 C 点和 C 的顺时针下一个点 C+1 的方向相同。
射线的交点就是我们的 \(I_{b,c}\)。这里的交点怎么求取呢？
- 平行的时候，\(\vec{B+1, B}\) 与 \(\vec{C+1,C}\) 叉乘为 0.我们将交点 \(I_{b,c}\) 与 BC 之间的距离定为无穷大。
- 根据平行叉乘为 0，以及克莱默公式有：
- 通过面积最大值来得到 3-stable 三角形。
于是我们有：
auto RotateAndKill = [&](int a, int b, int c) -> double {int ans_a = 1, ans_b = 2, ans_c = 3;double max_a = 0.0;int r = a, t = c;// Section 4:while (b != t || c != r) {// repeat 1. find A.while (dist(b, c, a + 1) > dist(b, c, a)) {a = a % stack_pt + 1;}double new_area = area(a, b, c);// output the 3-stable.if (new_area > max_a) {max_a = new_area;ans_a = a, ans_b = b, ans_c = c;}// Repeat 2: Calculate Ib,c.double dx1 = points[stack[b]][0] - points[stack[b + 1]][0], dy1 = points[stack[b + 1]][1] - points[stack[b]][1];double dx2 = points[stack[c]][0] - points[stack[c + 1]][0], dy2 = points[stack[c + 1]][1] - points[stack[c]][1];double delta = dy1 * dx2 - dy2 * dx1;if (delta >= -1e-9)b = b % stack_pt + 1;else {double U = dy1 * points[stack[c]][0] + dx1 * points[stack[c]][1];double V = dy2 * points[stack[b]][0] + dx2 * points[stack[b]][1];double Ix = (U * dx2 - V * dx1) / delta;double Iy = (V * dy1 - U * dy2) / delta;// (yc - yb) * xa + (xb - xc) * ya;if (dist(b, c, a) > (points[stack[c]][1] - points[stack[b]][1]) * Ix + (points[stack[b]][0] - points[stack[c]][0]) * Iy) {c = c % stack_pt + 1;} else {b = b % stack_pt + 1;}}}a = ans_a, b = ans_b, c = ans_c;return max_a;};
真的很有意思。花了我接近 6 个小时。