首先这是要解决什么问题:一个带权完全图,每条边都有自己的花费值cost[i]和收益值benifit[i],如果用x[i]来代表一条边取或不取,那么求一个生成树。要求:r=(∑cost[i]*x[i] ) / (∑benifit[i]*x[i] )最小。
经典题目:POJ2728 - Desert King
如何来求解:这里用到了0-1分数规划思想,对于上式可以变形为 z(r)=∑cost[i]*x[i] -r*∑benifit[i]*x[i]。而z(r)=0为我们所求。
这里有个非常重要的结论:z(r)为单调递减函数,因此是线性的。于是"我们可以兴高采烈地把z(r)看做以 cost[i]-r*benifit[i] 为边权的最小生成树的总权值"。显然,我们所求即为z(max(r))=0。这里max(r)是由z=0确定的。
于是有了两种算法:
1、Dinkelbanch算法,即迭代法。取r的一个初值,一般为0,得到一个生成树之后,令r=(∑benifit[i]*x[i] ) / (∑cost[i]*x[i] ),哐哐迭代之后……r是趋于最大值的,而z趋于0。
2、二分法,z(r)=0是我们所求,而它又是单调函数,那么二分r即可。据说二分比迭代慢了很多。
模板Dinkelbach POJ 2728:
#include <iostream> #include <cstdio> #include <cmath> #include <vector> #include <cstring> #include <algorithm> #include <string> #include <set> #include <ctime> #include <queue> #include <map> #include <sstream> #define CL(arr, val) memset(arr, val, sizeof(arr)) #define REP(i, n) for((i) = 0; (i) < (n); ++(i)) #define FOR(i, l, h) for((i) = (l); (i) <= (h); ++(i)) #define FORD(i, h, l) for((i) = (h); (i) >= (l); --(i)) #define L(x) (x) << 1 #define R(x) (x) << 1 | 1 #define MID(l, r) (l + r) >> 1 #define Min(x, y) x < y ? x : y #define Max(x, y) x < y ? y : x #define E(x) (1 << (x)) const double eps = 1e-6; typedef long long LL; using namespace std; const int N = 1<<10; const int inf = ~0u>>2; struct node { int x, y, z; } point[N]; double cost[N][N]; double benifit[N][N]; double mi[N]; int pre[N]; bool vis[N]; double c, d; int n; double prim(double ans) { int i, j, f; double mx; c = d = 0; for(i = 1; i < n; ++i) { mi[i] = cost[0][i] - ans*benifit[0][i]; pre[i] = 0; vis[i] = false; } mi[0] = 0; vis[0] = true; for(i = 0; i < n - 1; ++i) { mx = inf; f = 0; for(j = 0; j < n; ++j) { if(!vis[j] && mx > mi[j]) { mx = mi[j]; f = j; } } c += cost[pre[f]][f]; d += benifit[pre[f]][f]; vis[f] = true; for(j = 1; j < n; ++j) { double val = cost[f][j] - ans*benifit[f][j]; if(!vis[j] && mi[j] > val) { pre[j] = f; mi[j] = val; } } } return c/d; } int iabs(int x) { return x < 0 ? -x : x; } void init() { int i, j; for(i = 0; i < n; ++i) { for(j = 0; j < n; ++j) { cost[i][j] = iabs(point[i].z - point[j].z); benifit[i][j] = sqrt(1.*(point[i].x - point[j].x)*(point[i].x - point[j].x) + 1.*(point[i].y - point[j].y)*(point[i].y - point[j].y)); } } } void solve() { double ans = 0, tmp; while(1) { tmp = prim(ans); if(fabs(ans - tmp) < eps) break; else ans = tmp; } printf("%.3f\n", ans); } int main() { freopen("data.in", "r", stdin); int i; while(scanf("%d", &n), n) { for(i = 0; i < n; ++i) { scanf("%d%d%d", &point[i].x, &point[i].y, &point[i].z); } init(); solve(); } return 0; }