内接于半圆的最大矩形

在平面几何中，半圆常常出现在各种问题中。而求解内接于半圆的最大矩形也是一个经典的问题。

问题描述

给定一个半圆的圆心坐标 $(x_0, y_0)$ 和半径 $r$，求内接于该半圆的最大矩形。最大矩形需要满足其中两个顶点坐标在半圆上，其余两个顶点坐标在半圆内部。

解决方案

思路一：暴力枚举

暴力枚举是最朴素的解决方案。对于每个矩形，判断其是否满足题目要求，并记录面积最大的矩形。

时间复杂度为 $O(n^4)$，无法通过本题。

思路二：枚举对角线

观察矩形，发现其可以使用两条对角线确定，因此可以枚举所有对角线，然后判断对应的矩形是否符合条件。

在判断矩形是否符合条件时，可以计算出另外两个顶点的坐标并判断其是否在半圆内部。

时间复杂度为 $O(n^3)$，可以通过本题。

思路三：三分+计算几何+动态规划

由于所有顶点都在半圆上，可以通过枚举矩形的一条边，然后在半圆上二分查找与其垂直的另一条边的位置，进而确定这个矩形。

具体而言，可以先枚举矩形的一条边 $AB$，这条边可以在半圆上沿着圆心角来进行枚举。对于每个枚举的角度 $theta$，可以计算出直线 $l$ 的斜率 $k$，然后将 $l$ 与圆心 $O$ 的连线平移一段距离，就可以得到垂线 $l'$。接下来可以用三分查找来计算出垂线 $l'$ 与半圆的交点，从而确定点 $B'$。

最后可以使用动态规划来求解最大矩形。设 $dp_{i,j}$ 表示以点 $i$ 作为矩形的一条边，点 $j$ 作为矩形的对角线另一点，所能达到的最大面积。分类讨论如下：

• 如果 $i$ 和 $j$ 无法组成合法的矩形，则 $dp_{i,j}=0$；
• 如果 $i$ 和 $j$ 能组成合法的矩形，则对于所有 $k \in [i+1, j-1]$，有 $dp_{i,j}=\max{dp_{i,k}+dp_{k,j}+s_{i,k,j}}$，其中 $s_{i,k,j}$ 表示以 $i,j$ 和 $k$ 三个点构成的小三角形的面积。

最终的答案即为 $\max{dp_{i,j}}$。

时间复杂度为 $O(n^2\log n)$，可以通过本题。

代码实现

思路二实现

struct Point {
    double x, y;
    Point() {}
    Point(double x, double y): x(x), y(y) {}
};

struct Circle {
    Point o;
    double r;
    Circle() {}
    Circle(Point o, double r): o(o), r(r) {}
};

const double pi = acos(-1.0);
const double eps = 1e-8;

inline int sgn(double x) {
    return fabs(x) < eps ? 0 : (x < 0 ? -1 : 1);
}

inline double dist(Point p1, Point p2) {
    return hypot(p1.x - p2.x, p1.y - p2.y);
}

// 计算垂线交半圆的两个交点
inline vector<Point> get_foot_of_perpendicular(Circle C, Point P, double theta) {
    vector<Point> res;
    double r = C.r, x0 = C.o.x, y0 = C.o.y, x1 = P.x, y1 = P.y;

    double k = tan(theta);
    double b = y1 - k * x1;

    double A = k * k + 1;
    double B = -2 * k * b - 2 * x0 + 2 * y0 * k;
    double C = x0 * x0 + (b - y0) * (b - y0) - r * r;

    double delta = B * B - 4 * A * C;

    if (sgn(delta) < 0) {
        return res;
    }

    double x2 = (-B + sqrt(delta)) / (2 * A);
    double x3 = (-B - sqrt(delta)) / (2 * A);

    double y2 = k * x2 + b;
    double y3 = k * x3 + b;

    res.push_back(Point(x2, y2));
    res.push_back(Point(x3, y3));

    return res;
}

// 判断点是否在半圆内
inline bool in_circle(Circle C, Point P) {
    double r = C.r, x0 = C.o.x, y0 = C.o.y, x1 = P.x, y1 = P.y;
    return sgn(dist(C.o, P) - r) <= 0;
}

// 计算矩形面积
inline double calc_area(Point A, Point B, Point C) {
    double a = dist(B, C);
    double b = dist(A, C);
    return a * b;
}

double solve(Circle C) {
    double r = C.r, x0 = C.o.x, y0 = C.o.y;

    double ans = 0;

    // 枚举矩形的一条边
    for (double t1 = 0; t1 <= pi; t1 += pi / 200) {
        Point A(x0 + r * cos(t1), y0 + r * sin(t1));

        // 枚举矩形的对角线
        for (double t2 = t1; t2 <= pi; t2 += pi / 200) {
            Point B(x0 + r * cos(t2), y0 + r * sin(t2));

            // 计算垂线
            double k = tan(pi / 2 - t2);
            Point C(x0, y0 - 1);
            C.x = (B.y - A.y + k * A.x - k * B.x) / (2 * k);
            C.y = k * C.x - k * A.x + A.y;

            // 判断垂线交圆的两个交点是否在半圆内
            auto foots = get_foot_of_perpendicular(Circle(C, r), A, t1);
            if (!foots.empty() && in_circle(Circle(C, r), foots[0])) {
                Point D = foots[0];

                double area = calc_area(A, B, D);
                ans = max(ans, area);
            }

            foots = get_foot_of_perpendicular(Circle(C, r), B, t1);
            if (!foots.empty() && in_circle(Circle(C, r), foots[0])) {
                Point D = foots[0];

                double area = calc_area(A, B, D);
                ans = max(ans, area);
            }
        }
    }

    return ans;
}

思路三实现

struct Point {
    double x, y;
    Point() {}
    Point(double x, double y): x(x), y(y) {}
};

const double pi = acos(-1.0);
const double eps = 1e-8;

inline int sgn(double x) {
    return fabs(x) < eps ? 0 : (x < 0 ? -1 : 1);
}

inline double dist(Point p1, Point p2) {
    return hypot(p1.x - p2.x, p1.y - p2.y);
}

struct Circle {
    Point o;
    double r;
    Circle() {}
    Circle(Point o, double r): o(o), r(r) {}
};

// 计算两圆交点
inline vector<Point> get_intersection_points(Circle C1, Circle C2) {
    vector<Point> res;
    double r1 = C1.r, x1 = C1.o.x, y1 = C1.o.y, r2 = C2.r, x2 = C2.o.x, y2 = C2.o.y;

    double d = dist(C1.o, C2.o);

    if (sgn(d - r1 - r2) > 0 || sgn(d - fabs(r1 - r2)) < 0) {
        return res;
    }

    double a = (r1 * r1 - r2 * r2 + d * d) / (2 * d);
    double h = sqrt(r1 * r1 - a * a);
    double x0 = x1 + a * (x2 - x1) / d;
    double y0 = y1 + a * (y2 - y1) / d;

    res.push_back(Point(x0 + h * (y2 - y1) / d, y0 - h * (x2 - x1) / d));
    res.push_back(Point(x0 - h * (y2 - y1) / d, y0 + h * (x2 - x1) / d));

    return res;
}

// 计算向量 OA 和向量 OB 的叉积（OA 叉乘 OB）
inline double cross(Point O, Point A, Point B) {
    return (A.x - O.x) * (B.y - O.y) - (A.y - O.y) * (B.x - O.x);
}

// 判断点 C 是否在直线 AB 的左边（AB 向量逆时针旋转 90 度）
inline bool on_left(Point A, Point B, Point C) {
    return sgn(cross(A, B, C)) > 0;
}

// 计算三点组成的小三角形的面积（前提是 ABC 三点按逆时针顺序排列）
inline double calc_triangle_area(Point A, Point B, Point C) {
    return fabs(cross(A, B, C)) / 2;
}

// 计算矩形面积
inline double calc_area(Point A, Point B, Point C) {
    double a = dist(B, C);
    double b = dist(A, C);
    return a * b;
}

// 判断点是否在半圆内
inline bool in_circle(Circle C, Point P) {
    double r = C.r, x0 = C.o.x, y0 = C.o.y, x1 = P.x, y1 = P.y;
    return sgn(dist(C.o, P) - r) <= 0;
}

double solve(Circle C) {
    int n = 200;
    double r = C.r, x0 = C.o.x, y0 = C.o.y;

    double ans = 0;

    // 枚举矩形的一条边
    for (double t1 = 0; t1 <= pi; t1 += pi / n) {
        Point A(x0 + r * cos(t1), y0 + r * sin(t1));

        // 三分查找矩形的对角线
        for (double t2 = t1 + pi / n; t2 <= pi; t2 += pi / n) {
            Point B(x0 + r * cos(t2), y0 + r * sin(t2));

            double l = t1 + pi / 2, r = t2 + pi / 2;

            for (int _ = 0; _ < 50; ++_) {
                double m1 = l + (r - l) / 3;
                double m2 = l + (r - l) / 3 * 2;

                Point P1(x0 + r * cos(m1), y0 + r * sin(m1));
                Point P2(x0 + r * cos(m2), y0 + r * sin(m2));

                Circle C1(Circle(C.o, r), P1);
                Circle C2(Circle(C.o, r), P2);

                auto ps1 = get_intersection_points(Circle(C.o, r), C1);
                auto ps2 = get_intersection_points(Circle(C.o, r), C2);

                if (ps1.empty() || ps2.empty()) {
                    continue;
                }

                Point D1 = ps1.front(), D2 = ps2.front();
                if (in_circle(Circle(C.o, r), D1) && in_circle(Circle(C.o, r), D2)) {
                    double area1 = calc_triangle_area(A, B, D1);
                    double area2 = calc_triangle_area(A, B, D2);
                    double area3 = calc_area(A, B, D1);
                    double area4 = calc_area(A, B, D2);

                    ans = max(ans, area1 + area2 + area3 + area4);
                    l = m1;
                } else {
                    r = m2;
                }
            }
        }
    }

    return ans;
}