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;
}
//8092K	188MS    Dinkelbach  algorithm
//8092K	1125MS   二分