poj2728Desert King(最优比率生成树)

题目请戳这里

题目大意：一个王国n个城市编号1-n，1是首都，现在每个城市有坐标(x,y)和海拔h，现在在n个城市直接修水渠，要保证每个城市有水，水渠是水平的，但是由于每个城市的海拔不同，所以2个城市间修水渠的时候要一个cost去修水泵，cost=要修水渠的城市之间的海拔高度差，现在要求修单位长度的水渠的cost最小。

题目分析：n个城市修水渠，最少要修n-1条水渠，就是求给定图的一个生成树。同时生成树要保证sigma(cost)/sigma(len)最小。

这是一个最优比率生成树的问题。

详细论文戳这里

最优比率生成树的解法一般2种：二分法和Dinkelbach算法。

二分法在上面的论文里讲的非常详细了。构造的函数满足单调递减的，并且在最值处为0，二分还是比较容易的，但是二分的时间效率不高，所以这里稍微说一下Dinkelbach算法：

Dinkelbach算法其实就是一个牛顿迭代法，牛顿迭代的收敛速度是比较快的，也是Dinkelbach算法高效的原因。

上文中构造了一个子问题：

z(l) = cx' - l * dx'，其中l是迭代值，x'是一个子问题z(l)的最优解，因为每次都是求的最小生成树，那么观察z(l)这条直线，这条直线与z(L)在L= l处相切，并且与L轴相交于cx'/dx'，那么将cx'/dx'作为下一个迭代值，继续迭代下去，就会不断逼近最大值L*，因为上面论文已经证明z(L*) = 0。

关于Dinkelbach算法的详细分析请戳这里和这里

本题由于任意2点都是连通的，所以是一个及其稠密的图，所以在求mst的时候注意用prim，第一次写prim，用优先队列优化的，竟然TLE了。。。

直接暴力更新的反而跑得快。。。

详情请见代码：

#include <iostream>#include<cstdio>#include<cstring>#include<algorithm>#include<cmath>#include<queue>using namespace std;const int N = 1005;const int M = 1000005;const double INF = 1000000000.0;const double eps = 1e-5;struct node{    int x,y,h;}vex[N];bool vis[N];int last[N];double dis[N][N];double lowcost[N];int n;double getdis(struct node a,struct node b){    return sqrt((double)(a.x - b.x) * (a.x - b.x) + (double)(a.y - b.y) * (a.y - b.y));}double prim(double x){    memset(vis,false,sizeof(vis));    double retdis = 0;    int retcost = 0;    int sel = n - 1;    vis[1] = true;    int i;    for(i = 2;i <= n;i ++)    {        lowcost[i] = abs(vex[1].h - vex[i].h) - x * dis[1][i];        last[i] = 1;    }    while(sel --)    {        double Min = INF;        int u;        for(i = 1;i <= n;i ++)        {            if(vis[i] == true)                continue;            if(lowcost[i] < Min)            {                Min = lowcost[i];                u = i;            }        }        vis[u] = true;        retcost += abs(vex[last[u]].h - vex[u].h);        retdis += dis[last[u]][u];        for(i = 1;i <= n;i ++)        {            double cal = abs(vex[u].h - vex[i].h) - x * dis[u][i];            if(vis[i] == false && lowcost[i] > cal)            {                lowcost[i] = cal;                last[i] = u;            }        }    }    //return retcost - x * retdis;二分返回值    return (double)retcost/retdis;}int main(){    int i,j;    while(scanf("%d",&n),n)    {        for(i = 1;i <= n;i ++)        {            scanf("%d%d%d",&vex[i].x,&vex[i].y,&vex[i].h);        }        for(i = 1;i <= n;i ++)        {            for(j = 1;j <= i;j ++)            {                if(i == j)                    dis[i][j] = 0;                else                {                    dis[i][j] = dis[j][i] = getdis(vex[i],vex[j]);                }            }        }        double tmp,ans = 0;//        double l,r,mid;二分代码，注意精度//        l = 0.0;r = 100.0;//        while(r - l > eps)//        {//            mid = (l + r)/2;//            tmp = prim(mid);//            if(tmp < eps)//                r = mid;//            else//                l = mid;//        }//        ans = r;        while(1)        {            tmp = prim(ans);            if(fabs(ans - tmp) < eps)                break;            ans = tmp;        }        printf("%.3f\n",ans);    }    return 0;}//8092K188MS    Dinkelbach  algorithm//8092K1125MS   二分