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=(∑cost)÷(∑dist)

所以对于其他选择，一定有

(∑costi)÷(∑disti)≥rmin

稍微变形一下：

(∑costi)−(∑disti)∗rmin≥0(1)

因为dist一定是正数，所以我们构造的这么一个函数具有单调递减的性质：

rate（r）=(∑costi)−(∑disti)∗r(2)

根据此题的题意，我们需要找到使得rate函数为0的那个r

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

然后注意，我们在选定了r的数值之后，可以通过r，构造出新的边权：costi−disti∗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=(∑costi−rate(r))／(∑disti)(3)

构造一个新的 r2 令

r2=(∑costi)／(∑disti)(4)

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

当rate(r) < 0 时，r2<r,我们上一次构造出来的r不够小，所以下一次用r2

重复以上迭代的过程，最终能找到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);
}
}