poj 2728 Desert King （最优比率生成树/01分数规划）

思路：

更新于2017.7.6：
思路二的迭代算法貌似被称为：Dinkelbach
然后有关于这两种方法的选用，速度差距并不是特别大，我做过的这几道题里最多才差4倍（poj2728 二分1300ms dinkelbach 300ms limit 3000ms）
但是有时题目对dinkelbach极其不友好，因为我们很难得到suma和sumb，所以也就不容易构造出r。这时候二分的优势就体现出来了，因为它只需判断就好了。
poj3621就很好的体现了这点。（判断是否有负环容易，但要是取这个环的新权值就很麻烦了，因为在spfa过程中是不关心哪个点更新的距离，只关心更新的距离。。）
poj3621连接：
http://blog.csdn.net/wing_wuchen/article/details/74536630

====================================

思路1，01规划+二分：

设最小的比率为rmin，花费为cost，距离为dist，即有：

$$rmin = (\sum cost) \div (\sum dist)$$
所以对于其他选择，一定有
$$(\sum cost_i) \div (\sum dist_i) \ge rmin$$
稍微变形一下：
$$(\sum cost_i) - (\sum dist_i) * rmin \ge 0 \tag1$$
因为dist一定是正数，所以我们构造的这么一个函数具有单调递减的性质：
$$rate（r） = (\sum cost_i) - (\sum dist_i) * r \tag2$$
根据此题的题意，我们需要 找到使得rate函数为0的那个r

所以我们可以二分查找（0～100000000）（但其实r取到100就够了）

然后注意，我们在选定了r的数值之后，可以通过r，构造出新的边权：

$$cost_i - dist_i * r$$
在这个边权的基础上我们想要找到最小的，所以最小生成树prim上～（这个图是完全图，prim复杂度只与点的个数有关）

当然如果1式你构造的是反着的话，下边对应的求最大生成树就好。

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <math.h>
#include <queue>
#include <string.h>

using namespace std;

struct node{
    int x,y,z;
}a[1010];
int n;

double mapp[1010][1010];
double cost[1010][1010];
int vis[1010];
double dis[1010];

double discal(int i,int j){
    return sqrt( (a[i].x-a[j].x)*(a[i].x-a[j].x) + (a[i].y-a[j].y)*(a[i].y-a[j].y) );
}
double costcal(int i,int j){
    return abs(a[i].z-a[j].z);
}

const double inf = 1e15;
double rate(double r){
    memset(vis,0,sizeof(vis));
    for(int i = 1;i <= n;i++){
        dis[i] = inf;
    }
    double sum = 0;
    dis[1] = 0;
    for(int i = 1;i <= n;i++){
        int mark = -1;
        for(int j = 1;j <= n;j++){
            if(!vis[j]){
                if(mark == -1){
                    mark = j;
                }
                else if(dis[j] < dis[mark]){
                    mark = j;
                }
            }
        }
        if(mark == -1) break;
        vis[mark] = 1;
        sum += dis[mark];
        for(int j = 1;j <= n;j++){
            double temp = cost[mark][j] - mapp[mark][j] * r;
            if(!vis[j] && dis[j] > temp){
                dis[j] = temp;
            }
        }
    }
    return sum;
}

int main(){
    while(scanf("%d",&n),n){
        for(int i = 1;i <= n;i++){
            scanf("%d%d%d",&a[i].x,&a[i].y,&a[i].z);
        }
        for(int i = 1;i <= n;i++){
            for(int j = i+1;j <= n;j++){
                mapp[i][j] = mapp[j][i] = discal(i,j);
                cost[i][j] = cost[j][i] = costcal(i,j);
            }
        }
        double r = 100;
        double l = 0,mid;
        while(r - l > 1e-6){
            mid = (l+r)/2;
            double temp = rate(mid);
            if(temp > 0){
                l = mid;
            }
            else{
                r = mid;
            }
        }
        printf("%.3f\n",mid);
    }
}

思路二，01规划+迭代

将2式变形：

$$r = ( \sum cost_i - rate(r) )／(\sum dist_i) \tag3$$
构造一个新的 $r_2$ 令
$$r_2 = ( \sum cost_i )／(\sum dist_i) \tag4$$

然后我们发现当rate(r) > 0 时，$r_2 > r$,我们上一次构造出来的r不够大，所以下一次用$r_2$

当rate(r) < 0 时，$r_2 < r$,我们上一次构造出来的r不够小，所以下一次用$r_2$
重复以上迭代的过程，最终能找到1e-6误差范围内的结果。（即rate()这个函数收敛于rmin）

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <math.h>
#include <queue>
#include <string.h>

using namespace std;

struct node{
    int x,y,z;

}a[1010];
int n;

double mapp[1010][1010];
double cost[1010][1010];
int vis[1010];
double dis[1010];
int pre[1010];

double discal(int i,int j){
    return sqrt( (a[i].x-a[j].x)*(a[i].x-a[j].x) + (a[i].y-a[j].y)*(a[i].y-a[j].y) );
}
double costcal(int i,int j){
    return abs(a[i].z-a[j].z);
}

const double inf = 1e15;
double rate(double r){
    memset(vis,0,sizeof(vis));
    for(int i = 1;i <= n;i++){
        dis[i] = inf;
        pre[i] = 0;
    }
    double sum = 0,sumcost = 0;
    dis[1] = 0;
    for(int i = 1;i <= n;i++){
        int mark = -1,cos;
        for(int j = 1;j <= n;j++){
            if(!vis[j]){
                if(mark == -1){
                    mark = j;
                }
                else if(dis[j] < dis[mark]){
                    mark = j;
                }
            }
        }
        if(mark == -1) break;
        vis[mark] = 1;
        sum += mapp[pre[mark]][mark];//zheliguile
        //sum += dis[mark] //dis里存放的是构造出来的新的边权的值 
        sumcost += cost[pre[mark]][mark];
        for(int j = 1;j <= n;j++){
            double temp = cost[mark][j] - mapp[mark][j] * r;
            if(!vis[j] && dis[j] > temp){
                dis[j] = temp,pre[j] = mark;
            }
        }
    }
    return sumcost / sum;//return sum 返回的是新的r
}

int main(){
    while(scanf("%d",&n),n){
        for(int i = 1;i <= n;i++){
            scanf("%d%d%d",&a[i].x,&a[i].y,&a[i].z);
        }
        for(int i = 1;i <= n;i++){
            for(int j = i+1;j <= n;j++){
                mapp[i][j] = mapp[j][i] = discal(i,j);
                cost[i][j] = cost[j][i] = costcal(i,j);
            }
        }
        double r = 0;
        while(fabs(r-rate(r)) >  1e-6){
            r = rate(r);
        }
        printf("%.3f\n",r);
    }
}