QOJ4886 Best Sun 题解

First Post:

Last Update:

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

这题显然是要二分答案的，先二分一个 ，然后最大化 。

考虑如何 DP 中间那一层的凸包，并计算外面那些点的贡献。首先要枚举凸包 坐标最小 的点 ，然后将所有点对 按照向量 进行极角排序。

这里要讲极角排序的正确姿势：先按象限排序，再按叉积排序。也就说，对向量 ，定义布尔函数 ，即其是否在 轴右半边，在比较向量 时，如果 ，那么 大的要排在前面，否则看 的值，如果 ，那么 就应该排在前面，反之亦然。

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

先不考虑外面点的贡献，那么枚举已经排好序的点对 。如果 围成的三角形中没有包含其他的 ， 直接进行转移： 接下来考虑外面点的贡献怎么加入：

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

如果其投影在某条线段 上，这一部分贡献是可以预处理的，直接枚举 即可。需要判断 是否在向量 的右边，以及用点积判断其在直线 上的投影是否在线段 上。

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

以点 为例，我们考虑向量 与 轴的夹角 ，以及 周围两条边与 轴的夹角 ，容易发现 是符合条件的点当且仅当 。即 。统计这一部分贡献前，我们将所有点对 按照 （ 表示将向量 逆时针旋转 的结果）再排一遍序。然后在扫描按 排序的数组时 ，再用一个指针扫描这个按 排序的数组，把扫到的点加入贡献即可。

于是一次 固定的 Check 就可以在 的时间内解决，只要我们能 查询一个 中是否存在其它 中的点 。

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
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;
}

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

这题终于做完了，时间复杂度为 。

未尽细节参见代码。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
#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
*/
≡