主要参考:http://blog.youkuaiyun.com/sdj222555/article/details/7490797
http://www.cnblogs.com/lotus3x/archive/2009/03/21/1418480.html
概念
有带权图G, 对于图中每条边e[i], 都有benifit[i](收入)和cost[i](花费), 我们要求的是一棵生成树T, 它使得 ∑(benifit[i]) / ∑(cost[i]), i∈T 最大(或最小).
这显然是一个具有现实意义的问题.
解法之一 0-1分数规划
设x[i]等于1或0, 表示边e[i]是否属于生成树.
则我们所求的比率 r = ∑(benifit[i] * x[i]) / ∑(cost[i] * x[i]), 0≤i<m .
为了使 r 最大, 设计一个子问题---> 让 z = ∑(benifit[i] * x[i]) - l * ∑(cost[i] * x[i]) = ∑(d[i] * x[i]) 最大 (d[i] = benifit[i] - l * cost[i]) , 并记为z(l). 我们可以兴高采烈地把z(l)看做以d为边权的最大生成树的总权值.
然后明确两个性质:
1. z单调递减
证明: 因为cost为正数, 所以z随l的减小而增大.
2. z( max(r) ) = 0
证明: 若z( max(r) ) < 0, ∑(benifit[i] * x[i]) - max(r) * ∑(cost[i] * x[i]) < 0, 可化为 max(r) < max(r). 矛盾;
若z( max(r) ) >= 0, 根据性质1, 当z = 0 时r最大.
其实原题就是求 MIN( ∑CiXi / ∑DiXi ) Xi∈{0,1} ,对每个生成树,设其比率r=∑CiXi / ∑DiXi ,可得∑CiXi - ∑DiXi * r=0(条件1)
那么对于所有的生成树,显然∑CiXi - ∑DiXi * min(r) >= 0,当 ∑CiXi / ∑DiXi = min(r)时,等号成立。 而我们现在不知道min(r)是多少,只好进行枚举,对每个枚举的r ,构建新的权值(Ci-Di*r),然后求最小生成树, 为什么求最小呢? 我的理解就是这是为了寻找使得生成树的总权值为0的可能性,因为只有当其等于0 的时候,才满足了条件1 这个条件, 说明这个r是可行的,并且如果r枚举到值为min(r)时,其最小生成树的的总权值必然恰好等于0,但是如果不能等于0, 比如大于0, 显然是对该r值,所有的生成树上无论如何也满足不了条件1,说明r值就是偏小了。同理如果小于0,r值是偏大的,说明可能存在某些生成树使得满足条件1,而我们的目标是在满足条件1的情况下使得r最小。
根据这个我们可以发现,实际上r的值是可以进行二分查找的。 而也有人给出了更为高效的迭代方法。
第一种是二分的方法
#include <iostream>
#include <cstdio>
#include <cstring>
#include <vector>
#include <stack>
#include <cmath>
#define INF 1000000000
#define eps 1e-7
using namespace std;
const int N = 1005;
int n;
double e[N][N],low[N];
int nearn[N];
struct point{
double x,y,z;
}p[N];
double getdis(point a, point b){
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
bool prim(double l){
double sum = 0;
for(int i = 1; i <= n; i++){
nearn[i] = 1;
low[i] = abs(p[1].z - p[i].z) - e[1][i] * l;
}
nearn[1] = -1;
for(int i = 1; i < n; i++){
double mi = INF;
int v = -1;
for(int j = 1; j <= n; j++)
if(nearn[j] != -1 && low[j] < mi){
v = j;
mi = low[j];
}
if(v == -1)continue;
nearn[v] = -1;
sum += low[v];
for(int j = 1; j <= n; j++){
double tmp = abs(p[v].z - p[j].z) - e[v][j] * l;
if(nearn[j] != -1 && tmp < low[j])
{
low[j] = tmp;
nearn[j] = v;
}
}
}
return sum >= 0;
}
int main() {
while(~scanf("%d", &n) && n) {
for(int i = 1; i <= n; i++)
scanf("%lf%lf%lf", &p[i].x, &p[i].y, &p[i].z);
for(int i = 1; i <= n; i++)
for(int j = 1; j <= n; j++)
e[i][j] = getdis(p[i], p[j]);
double l = 0.0, r = 50.0, mid;
while(r - l > eps) {
mid = (l + r) / 2;
if(prim(mid)) l = mid;
else r = mid;
}
printf("%.3f\n", r);
}
return 0;
}
第二种是迭代
#include <iostream>
#include <cstdio>
#include <cstring>
#include <vector>
#include <stack>
#include <cmath>
#define INF 1000000000
#define eps 1e-7
using namespace std;
const int N = 1005;
int n;
double e[N][N],low[N];
int nearn[N];
struct point{
double x,y,z;
}p[N];
double getdis(point a, point b){
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
double prim(double l){
double cost = 0, len = 0;
double sum = 0;
for(int i = 1; i <= n; i++){
nearn[i] = 1;
low[i] = abs(p[1].z - p[i].z) - e[1][i] * l;
}
nearn[1] = -1;
for(int i = 1; i < n; i++){
double mi = INF;
int v = -1;
for(int j = 1; j <= n; j++)
if(nearn[j] != -1 && low[j] < mi){
v = j;
mi = low[j];
}
if(v == -1)continue;
cost += abs(p[nearn[v]].z - p[v].z);
len += e[nearn[v]][v];
nearn[v] = -1;
sum += low[v];
for(int j = 1; j <= n; j++){
double tmp = abs(p[v].z - p[j].z) - e[v][j] * l;
if(nearn[j] != -1 && tmp < low[j])
{
low[j] = tmp;
nearn[j] = v;
}
}
}
return cost/len;
}
int main() {
while(~scanf("%d", &n) && n) {
for(int i = 1; i <= n; i++)
scanf("%lf%lf%lf", &p[i].x, &p[i].y, &p[i].z);
for(int i = 1; i <= n; i++)
for(int j = 1; j <= n; j++)
e[i][j] = getdis(p[i], p[j]);
double a = 0, b;
while(true){
b = prim(a);
if(fabs(a-b) < eps) break;
a = b;
}
printf("%.3f\n", a);
}
return 0;
}