qoj4886 Best Sun 题解

天命之路 · 2023-10-20 12:02:00 · 个人记录

本题细节极多，该题解包含了 DP 凸包的正确方式 以及解决另一个问题的正确方式。

这题显然是要二分答案的，先二分一个 mid，然后最大化 S - mid \times P。

考虑如何 DP 中间那一层的凸包，并计算外面那些点的贡献。首先要枚举凸包 x 坐标最小 的点 s ，然后将所有点对 (i,j) 按照向量 v_{i,j} = p_j - p_i 进行极角排序。

这里要讲极角排序的正确姿势：先按象限排序，再按叉积排序。也就说，对向量 v，定义布尔函数 sgn(v) = [v_x > 0 \lor (v_x = 0 \land v_y > 0)]，即其是否在 y 轴右半边，在比较向量 v_1,v_2 时，如果 sgn(v_1) \ne sgn(v_2)，那么 sgn 大的要排在前面，否则看 v_1 \times v_2 的值，如果 >0，那么 v_1 就应该排在前面，反之亦然。

我们以逆时针方向尝试填充这个凸包的每一个点，设 dp_i 表示当前凸包是 s \sim i 的一段凸壳时的最优答案（即当前最后一个填进去的点是 i）。那么最后只需检查 dp_s 是否 > 0。

先不考虑外面点的贡献，那么枚举已经排好序的点对 (i,j)。如果 p_s,p_i,p_j 围成的三角形中没有包含其他的 p_k，直接进行转移：

dp_j \gets \max\{dp_j,dp_i + S_{\triangle p_sp_ip_j} - mid \times |p_ip_j|\}

接下来考虑外面点的贡献怎么加入：

考虑一个外面点 p_k 及其在凸包上的投影。

如果其投影在某条线段 p_ip_j 上，这一部分贡献是可以预处理的，直接枚举 (i,j,k) 即可。需要判断 p_k 是否在向量 \vec{p_ip_j} 的右边，以及用点积判断其在直线 p_ip_j 上的投影是否在线段 p_ip_j 上。

如果其投影在某个顶点 p_i 上，事情会有点麻烦，我们需要计数有多少个多边形外的点的投影正好在顶点 p_i 上。准确地说，这些点在 p_i 边上两条线段上的投影都在这两条线段外。如下图，标黄的点就是我们在这一部分要统计的点。

以点 G 为例，我们考虑向量 \vec{BG} 与 x 轴的夹角 \beta，以及 p_i 周围两条边与 x 轴的夹角 \alpha,\gamma，容易发现 G 是符合条件的点当且仅当 \alpha - \dfrac{\pi}{2} < \beta < \gamma - \dfrac{\pi}{2}。即 \alpha < \beta + \dfrac{\pi}{2} < \gamma。统计这一部分贡献前，我们将所有点对 (i,j) 按照 v'_{i,j} = \operatorname{Rot}(p_j - p_i) （Rot(v) 表示将向量 v 逆时针旋转 \dfrac{\pi}{2} 的结果）再排一遍序。然后在扫描按 v_{i,j} 排序的数组时，再用一个指针扫描这个按 v' 排序的数组，把扫到的点加入贡献即可。

于是一次 s 固定的 Check 就可以在 O(n^2) 的时间内解决，只要我们能 O(1) 查询一个 \triangle p_sp_ip_j 中是否存在其它 p 中的点。

三角形是怎么求面积的？S_{\triangle ABC} = |\dfrac{A \times B + B \times C + C \times A}{2}|。其中，A \times B 实际上就是 S_{\triangle OAB} 的 2 倍（也可能是 -2 倍）。上述三个叉积加起来，正好可以把 \triangle ABC 外的所有部分全部抵消掉。考虑拓展这种思路，预处理 cnt_{i,j} ，如果一个点 k 严格处于 \triangle Op_ip_j 内部，让 cnt_{i,j} 加 2。如果其处于边界上，则让 cnt_{i,j} 加 1。在算 \triangle p_sp_ip_j 中包含的点数时，按照叉积的正负性给结果对应的贡献。这样，将结果取绝对值后，这个三角形每个内部点会被算两次，每个边界点会被算一次。因为我们只需知道三角形内点数是否为 0，我们直接判断这个算出来的值是不是 0 即可。这部分代码如下：

inline int NeedAdd(const Vec &A,const Vec &B,const Vec &p) { // 求出一个点应该算几次
    if(!InTriangle(O,A,B,p)) return 0;
    if(!p.x && !p.y) return 0;

    if((O - p) * (A - p) == 0 
        || (O - p) * (B - p) == 0 || (A - p) * (B - p) == 0) return 1;
    return 2;
}
inline bool CheckTriangle(int i,int j,int k) {

    if(InSegment(p[i],p[j],p[k])) return VL[i][j] <= 1; // VL[i][j] 表示处于线段 p_ip_j 上的点数，每个点算 1 次
    if(InSegment(p[i],p[k],p[j])) return VL[i][k] <= 1; 
    if(InSegment(p[j],p[k],p[i])) return VL[j][k] <= 1;
    int tmp = TagO - (p[i] == O) - (p[j] == O) - (p[k] == O);
    if(tmp && InTriangle(p[i],p[j],p[k],O)) return false;
    if(VL[i][j] || VL[j][k] || VL[i][k]) return false;

    int num = 0;
    int v1 = V[i][j] - NeedAdd(p[i],p[j],p[k]); // V[i][j] 表示前文提到的 cnt[i][j]
    int v2 = V[j][k] - NeedAdd(p[j],p[k],p[i]);
    int v3 = V[k][i] - NeedAdd(p[k],p[i],p[j]);
    if(p[i] * p[j] > 0) num += v1; else if(p[i] * p[j] < 0) num -= v1;
    if(p[j] * p[k] > 0) num += v2; else if(p[j] * p[k] < 0) num -= v2;
    if(p[k] * p[i] > 0) num += v3; else if(p[k] * p[i] < 0) num -= v3;

    return num == 0;
}

如果按照上文的写法，我们可以得到一个总复杂度为 O(n^3 \log V) 的程序（对于一个 s 的 Check 是 n^2 的，要对所有 s 都 Check），可能过不去。不妨设 ans_s 表示以 s 为起点的最终答案，然后对每个 s 分开二分。这样仍没有优化，我们不妨将 s 序列随机打乱，这样 ans_s 的前缀最大值就期望为 \log n 个。遍历到一个 s 时，我们可以用一次 Check 判断 ans_s 是否大于前面所有的 ans_i ，如果大于再进行二分。这样的期望时间复杂度就是 O(n^2 \log n \log V)。

这题终于做完了，时间复杂度为 O(n^3 + n^2 \log n \log V)。

未尽细节参见代码。

#include <bits/stdc++.h>
using namespace std;
const int N = 3e2 + 5;
typedef long long i64;
typedef double db;
typedef pair<int,int> pii;
#define FI first
#define SE second
const db eps = 1e-9,PI = acos(-1.0);
mt19937 Rnd(time(0));
struct Vec {
    int x,y;
    Vec(){}
    Vec(const int _x,const int _y):x(_x),y(_y){}
    inline Vec operator + (const Vec &rhs) const { return Vec(x + rhs.x,y + rhs.y);}
    inline Vec operator - (const Vec &rhs) const { return Vec(x - rhs.x,y - rhs.y);}
    inline i64 operator * (const Vec &rhs) const { return 1ll * x * rhs.y - 1ll * y * rhs.x;}
    inline bool operator == (const Vec &rhs) const { return x == rhs.x && y == rhs.y;}
};
Vec O(0,0);
inline bool Cmp(const Vec &x,const Vec &y) { return (x.x != y.x) ? (x.x < y.x) : (x.y < y.y);}
inline i64 Dot(const Vec &A,const Vec &B) { return 1ll * A.x * B.x + 1ll * A.y * B.y;}
inline db Distance(const Vec &A,const Vec &B) { return sqrt(Dot(A - B,A - B));}
inline db Area(const Vec &A,const Vec &B,const Vec &C) { return 0.5 * abs(A * B + B * C + C * A);}
inline bool Equal(const db &A,const db &B) { return fabs(A - B) < eps;}
inline bool InSegment(const Vec &A,const Vec &B,const Vec &p) { return Equal(Distance(A,p) + Distance(p,B),Distance(A,B));}
inline bool InTriangle(const Vec &A,const Vec &B,const Vec &C,const Vec &p) {
    if((B - A) * (p - A) >= 0 && (C - B) * (p - B) >= 0 && (A - C) * (p - C) >= 0) return true;
    if((B - A) * (p - A) <= 0 && (C - B) * (p - B) <= 0 && (A - C) * (p - C) <= 0) return true;
    return false;
}
inline Vec Rotate(const Vec &A) { return Vec(-A.y,A.x);}
inline int sgn(const Vec &A) { return A.x > 0 || A.x == 0 && A.y > 0;}
Vec p[N];
int n;
int V[N][N],VL[N][N]; 
int TagO;
db sum[N][N];
inline int NeedAdd(const Vec &A,const Vec &B,const Vec &p) {
    if(!InTriangle(O,A,B,p)) return 0;
    if(!p.x && !p.y) return 0;

    if((O - p) * (A - p) == 0 
        || (O - p) * (B - p) == 0 || (A - p) * (B - p) == 0) return 1;
    return 2;
}
inline bool CheckTriangle(int i,int j,int k) {

    if(InSegment(p[i],p[j],p[k])) return VL[i][j] <= 1;
    if(InSegment(p[i],p[k],p[j])) return VL[i][k] <= 1;
    if(InSegment(p[j],p[k],p[i])) return VL[j][k] <= 1;
    int tmp = TagO - (p[i] == O) - (p[j] == O) - (p[k] == O);
    if(tmp && InTriangle(p[i],p[j],p[k],O)) return false;
    if(VL[i][j] || VL[j][k] || VL[i][k]) return false;

    int num = 0;
    int v1 = V[i][j] - NeedAdd(p[i],p[j],p[k]);
    int v2 = V[j][k] - NeedAdd(p[j],p[k],p[i]);
    int v3 = V[k][i] - NeedAdd(p[k],p[i],p[j]);
    if(p[i] * p[j] > 0) num += v1; else if(p[i] * p[j] < 0) num -= v1;
    if(p[j] * p[k] > 0) num += v2; else if(p[j] * p[k] < 0) num -= v2;
    if(p[k] * p[i] > 0) num += v3; else if(p[k] * p[i] < 0) num -= v3;

    return num == 0;
}

struct Line {
    int s,t,type;
    Vec v;
    Line(){}
    Line(const int _s,const int _t,const int _type,const Vec _v):
        s(_s),t(_t),type(_type),v(_v){}
    inline bool operator < (const Line &rhs) const { 
        // return atan2(v.y,v.x) < atan2(rhs.v.y,rhs.v.x);
        int t1 = sgn(this->v),t2 = sgn(rhs.v);
        if(t1 != t2) return t1 > t2;
        i64 val = this->v * rhs.v;
        // i64 val = atan2(rhs.v.y,rhs.v.x) - atan2(v.y,v.x);
        if(val) return val > 0;
        else return type < rhs.type;
    }
};
Line L[N * N * 2];
int tot;
inline void Init() {
    cin >> n;
    for(int i = 1;i <= n;i++) cin >> p[i].x >> p[i].y;
    TagO = 0;
    for(int i = 1;i <= n;i++) if(p[i] == O) TagO++;
    sort(p + 1,p + n + 1,Cmp);
    for(int i = 1;i <= n;i++) 
        for(int j = 1;j <= n;j++)
            if(i < j)
            for(int k = 1;k <= n;k++) {
                if(k == i || k == j) continue;  
                if(p[k] == O) continue;
                if(!InTriangle(p[i],O,p[j],p[k])) continue;
                V[i][j] += NeedAdd(p[i],p[j],p[k]);
                if(InSegment(p[i],p[j],p[k])) ++VL[i][j];
            }
            else V[i][j] = V[j][i],VL[i][j] = VL[j][i];
    for(int i = 1;i <= n;i++)
        for(int j = 1;j <= n;j++)
            for(int k = 1;k <= n;k++) {
                if(k == i || k == j) continue;
                if((p[k] - p[i]) * (p[j] - p[i]) <= 0) continue;
                if(Dot(p[j] - p[i],p[k] - p[i]) >= 0 && Dot(p[i] - p[j],p[k] - p[j]) > 0)
                    sum[i][j] += min(Distance(p[k],p[i]),Distance(p[k],p[j]));
            }
    for(int i = 1;i <= n;i++)
        for(int j = 1;j <= n;j++) {
            if(i == j) continue;
            sum[i][j] += Distance(p[i],p[j]);
            L[++tot] = Line(i,j,0,p[j] - p[i]);
            L[++tot] = Line(i,j,1,Rotate(p[j] - p[i]));
        }
    sort(L + 1,L + tot + 1);
    cerr << 1.0 * clock() / CLOCKS_PER_SEC << endl;
}

db dp[N];
inline bool Check(int st,db mid) {
    for(int i = 1;i <= n;i++) dp[i] = -1e15;
    dp[st] = -eps;

    for(int i = 1;i <= tot;i++) {
        int s = L[i].s,t = L[i].t,tp = L[i].type;
        if(tp == 1) {
            if(s >= st) dp[s] -= mid * Distance(p[s],p[t]);
        } else 
            if(s >= st && t >= st && CheckTriangle(st,s,t))
                dp[t] = max(dp[t],dp[s] + Area(p[st],p[s],p[t]) - mid * sum[s][t]); 
    }

    return dp[st] > eps;
}
inline void Clr(int n) {
    for(int i = 1;i <= n;i++) for(int j = 1;j <= n;j++) V[i][j] = VL[i][j] = sum[i][j] = 0;
    tot = 0;
}
inline void work() {
    Clr(n);
    Init();
    static int perm[N];
    for(int i = 1;i <= n;i++) perm[i] = i;
    shuffle(perm + 1,perm + n + 1,Rnd);
    db ans = 0;

    for(int i = 1;i <= n;i++) 
        if(Check(perm[i],ans)) {
            db lef = ans,rig = 1e6;
            for(int _ = 1;_ <= 45;_++) {
                db mid = (lef + rig) / 2;
                if(Check(perm[i],mid)) lef = mid;
                else rig = mid;
            }
            ans = lef;
        }
    printf("%.10lf\n",ans);
    cerr << 1.0 * clock() / CLOCKS_PER_SEC << endl;
}
int main() {
    int T;
    cin >> T;
    while(T--) work();
    return 0;
}
/*
1
4
1 -4
3 -3
1 -1
2 1
*/