关于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: pop
            while (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: pop
            while (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 个小时。