前言
今天老师讲二分、迭代,顺便搜到了,结果刷了我一天...
题面
http://poj.org/problem?id=2728
题意
给定你N个点的坐标(x, y, z),对于任意两点Vi, Vj之间的连边,其长度为√[(xi - xj)2 + (yi - yj)2],花费为|zi - zj|。求这N个点的生成树所能达到的总花费和比总长度和的最小值。
分析
§ 1 大力猜结论
相信大家如果是头回看到这个题,一定会想去取花费比上长度的值尽可能小的边。但是这样是错误的。容易举出反例。
§ 2 转化问题
我们记这个最小值存在且为min
那么min = ∑Ci / ∑Li,这个由题意显然,其中Ci,Li分别为对应生成树上边的花费和长度。
我们考虑二分这个最小值,若当前找到一个值为x,则x要满足条件需要存在生成树使得
x ≥ ∑Ci / ∑Li
即 x∑Li ≥ ∑Ci
即 ∑(xLi - ∑Ci) ≥ 0
那么我们只需要每次找到一个x,然后把每条边的边权改为(xLi - ∑Ci),做一遍最大生成树,判断所需总价值是否大于等于0即可。
§ 3 变得更快
考虑这是一个完全图,那么当然用Prim做(我一开始写Kruskal了,too young & too simple...),考虑这个规模没有必要写堆优化(堆优化会超时,没错我又写过...)。
然后注意到Li和Ci是一直不会变的,所以不妨在二分前把他们预处理好,然后对于给定x,每条边的边权也可以在做生成树的时候再算出,无需预处理。这样会省掉很多冗余的运算,常数优秀。
§ 4 迭代
对呀!讲了这么多,还是在说二分嘛!
为什么用迭代是因为迭代比二分快很多。迭代,大概的原理就跟牛顿法过程差不多。
考虑我们二分,每次我们求出一个mid,检查一下答案,我们只判断了这个mid是否可行,然后用它来更新二分下界或上界。
但是对于可行解,打个比方,对于这个问题,我们求出的最大生成树其本身的这个总花费和比总长度和值可能就是比mid要优秀的,而我们接下来再二分时,却只利用了这个mid,却没有利用这个更加优秀的值。
这样就放弃了本来可以变得更小的答案区间。
迭代就是利用这个更加优秀的值(在解可行时),去更新这个答案区间。
想想牛顿迭代,这个速度必然是极其优秀的了!
§ 5 犯下愚蠢的错误
对于POJ这种只返回一个WA或者其他信息的清新简洁无广告良心网站,我们一定要学会利用好Discuss这个东西。
比如说这道题,幸好看了Discuss,要不然我估计今天都一直WA。
我用了double类型存放ans
用G++提交的情况下,用%.3lf,是WA的;而改成%.3f,A了!
其实double类型输入是用%lf,而输出是用%f的,我之间一直以为输出也是%lf...(不过好像据说是编译器的问题,G++里double就必须用%f输出)
参考程序
1 // POJ2728
2 #include <cstdio>
3 #include <cmath>
4 #include <cstring>
5 #include <algorithm>
6 const int MAXN = 1005;
7 const double EPS = 0.0001; // Discuss里说这个EPS要开小一点,不过这么多其实够了
8 namespace FastIO { // 快读
9 template <typename T>
10 inline bool read(T & x) {
11 x = 0;
12 char ch = getchar(); bool f = false, c = false;
13 for (; ch < '0' || ch > '9'; f |= (ch == '-'), ch = getchar());
14 for (; ch >= '0' && ch <= '9'; x = (x << 3) + (x << 1) + (ch ^ '0'), c = true, ch = getchar());
15 if (f) x = -x;
16 return c;
17 }
18
19 template <typename T>
20 inline void write(T & x) {
21 if (x == 0) return (void)(putchar('0'));
22 if (x < 0) { putchar('-'); x = -x; }
23 int num[15], len = 0;
24 for (; x; x /= 10) num[len++] = x % 10;
25 while (len) putchar(num[--len] ^ '0');
26 }
27
28 template <typename T>
29 inline void writeln(T & x) { write(x); putchar('
'); }
30 }
31
32 int N, E, X[MAXN], Y[MAXN], Z[MAXN], C[MAXN][MAXN], prv[MAXN];
33 double lb, ub, D[MAXN], L[MAXN][MAXN];
34 bool used[MAXN];
35
36 inline double sq(double x) { return x * x; }
37 void solve();
38 double check(double x);
39
40 int main() {
41 while (FastIO::read(N)) {
42 if (!N) break;
43 solve();
44 }
45 return 0;
46 }
47
48 void solve() {
49 using namespace FastIO;
50 int i, j;
51 for (i = 1; i <= N; i++) { read(X[i]); read(Y[i]); read(Z[i]); }
52 for (i = 1; i < N; i++) // 预处理一下长度和花费
53 for (j = i + 1; j <= N; j++) {
54 L[j][i] = L[i][j] = sqrt(sq(X[i] - X[j]) + sq(Y[i] - Y[j]));
55 C[j][i] = C[i][j] = std::abs(Z[i] - Z[j]);
56 ub = std::max(ub, C[i][j] / L[i][j]); // 找一个上界
57 }
58 double mid, tmp;
59 for (i = 0, lb = 0, ub += 1; i < 50 && ub - lb > EPS; i++) { // 迭代深度开50够了,其实100也不会TLE的
60 mid = (lb + ub) / 2;
61 tmp = check(mid);
62 if (tmp > 0) ub = tmp; // 如果解可行则用更优秀的值更新上界
63 else lb = mid; // 否则还是只能用mid更新下界
64 }
65 printf("%.3f
", round(ub * 1000) / 1000);
66 }
67
68 double check(double x) { // 裸的O(n^2)的Prim没什么好说的
69 int i, v;
70 double sum = 0, suml = 0, sumc = 0, tmp;
71 memset(used, 0, sizeof(used));
72 D[1] = 0; used[1] = true;
73 for (i = 2; i <= N; i++) D[i] = x * L[1][i] - C[1][i], prv[i] = 1; // 这个prv表示更新当前点距离的一个前驱,是为了找到那条边求出更优秀的值
74 while (true) {
75 v = -1;
76 for (i = 1; i <= N; i++)
77 if (!used[i] && (v == -1 || D[i] > D[v])) v = i;
78 if (v == -1) break;
79 used[v] = true; sum += D[v]; suml += L[prv[v]][v]; sumc += C[prv[v]][v];
80 for (i = 1; i <= N; i++)
81 if (!used[i]) {
82 tmp = x * L[v][i] - C[v][i];
83 if (tmp > D[i]) { D[i] = tmp; prv[i] = v; }
84 }
85 }
86 if (sum >= 0) return sumc / suml;
87 return -1;
88 }