• 网络又断了,我本来想再切两道题目的,这回安逸了,老天爷要给我放假,俄!天也!
  • frkstycArXoR两位大神侃的相当开心哈,过料,才发现这两个仁兄!
  • 看着标题以为很难,其实一点都不难,只要理论过关,两种方法都好过!
  • 理论参见:最优比率生成树分析[转一个兄弟的],关键是这个理论很重要,Dinkelbach算法。可以是说0-1分数规划的最强方法,很像牛顿哥哥的切线逼近法。

二分理论:z(rate)函数单减

1.z(r)<0 <=> r<r*

2.z(r)=0 <=> r=r*

3.z(r)>0 <=> r>r*

  • 下面用两种方法解两道题目,第一道是USACO US open 2001的2题 Earthquake

分析:题目要求r=(F-sum(CiXi))/(TiXi)的最大生成树,转化下就是求z(l)=F-sum(CiXi)-l*sum(TiXi)的最大值,然后再转化为z(l)=F-sum(CiXi+l*Sum(TiXi)).然后就是按照二分理论求解,关键在于这道题目是一个基于边的生成树,所以选择Kruskal算法。代码如下:

#include <iostream>
#include <math.h>
#include <fstream>

//UASCO earthquake
//1 <= N <= 400, 1 <= M <= 10,000, 1 <= F <= 2,000,000,000
//1 <= t <= 2,000,000,000, 1 <= c <= 2,000,000,000
using namespace std;
const int MAXN=400, MAXM=10000;
int n, m, G[MAXN+1][MAXN+1]={0};
long long f;
typedef struct Edge{
    int i, j; long double t, c, value;
}Edge, *Pedge;
Edge es[MAXM];
int eid=0, hsize;

//UFS implementation
class UFS{
private:
    int n[MAXN], setn;
public:
    UFS(){clr();}
    void clr(){memset((void*)n, 0, sizeof(n));}
    void setk(int k){
        this->setn=k; memset((void*)(n+1), 0xff, sizeof(int)*this->setn);
    }
    int find(int v){
        if(v<=0||v>setn) return -1;
        int p=v;
        while(n[v]>0) v=n[v];
        while(p!=v){
            int tmp=n[p];
            n[p]=v;
            p=tmp;
        }
        return v;
    }
    void unions(int i, int j){
        if(i<=0||i>setn||j<=0||j>setn) return;
        int v1=find(i), v2=find(j);
        if(v1==v2) return;
        int tmp1=n[v1], tmp2=n[v2];
        if(tmp1<tmp2){
            n[v1]=tmp1+tmp2;
            n[v2]=v1;
        }
        else{
            n[v2]=tmp1+tmp2;
            n[v1]=v2;
        }
    }
};
UFS _u;

//declaration
void MST(double l, long double& cc, long double& ct, bool inc);
//equation : di=l*ti+ci; 1. value=l*ti, l=1, 2. value=l*ti+ci.
int inline popHeap();

//implementation
void MST(double l, long double& cc, long double& ct, bool inc){//引用求出其和
    hsize=m; _u.setk(m);//for [1…m]
    for(int i=1; i<hsize+1; i++){//calculate the value of edge
        if(!inc) es[i].value=es[i].t*l;
        else es[i].value=es[i].t*l+es[i].c;
    }
    do{
        int tmp=popHeap(); int s=es[tmp].i, t=es[tmp].j;
        if(_u.find(s)!=_u.find(t)){
cc+=es[tmp].c, ct+=es[tmp].t, _u.unions(s, t);
        }
    }while(hsize>=1);
}
int inline popHeap(){
    if(hsize<1) return -1;
    int i=1;
    long double minValue=es[i].value;
    for(int j=2; j<hsize+1; j++){
        if(es[j].value<minValue) i=j, minValue=es[j].value;
    }
    if(i!=hsize){
        Edge tmp=es[i];es[i]=es[hsize];es[hsize]=tmp;
    }
    return hsize–;
}

#define fbase
int main()
{
#ifdef fbase
    ifstream in("input.txt");
    if(!in) return EXIT_FAILURE;
    cin.rdbuf(in.rdbuf());
#endif

    //implementation
    long double cc=0, ct=0, t=0, c=0, lr, hr, z, mid;
    double maxRate;
    cin>>n, cin>>m, cin>>f;
    for(int i=1; i<=m; i++){
        cin>>es[i].i, cin>>es[i].j, cin>>es[i].c, cin>>es[i].t;
        t+=es[i].t, c+=es[i].c;
    }
    lr=1e6*(double)(f-c)/(double)(t);//lower boundary
    MST(1,cc,ct,false);//with t as value
    hr=1e6*(double)f/(double)(ct);

    while(lr<hr){
        mid=lr+(hr-lr)/2;
        cc=0,ct=0,maxRate=mid*1e-6;//一个技巧问题,比率太小,精度应该适当扩大!
        MST(maxRate, cc, ct, true);
        z=f-(cc+maxRate*ct);
        if(fabs(z)<1e-6) break;//achieve the more exact value
        else if(z>0) lr=mid+1;
        else hr=mid-1;
    }
    if(maxRate<0) printf("0.0000");
    else printf("%.4f\n", maxRate);

#ifdef fbase
    in.close();
#endif
    return 0;
}
#undef fbase

  • 第二道题目,呵呵,楼天成哥哥扬名之作!当然,他肯定秒杀没问题!一个完全图Prim算法+Dinkebach算法。这个算法很巧妙,等于我们先用树生成算法得到一个估计值rate,然后用得到的那个向量来求下一个逼近值,直到rate切线重合![BOP的体现,下面是用数组实现的Prim+Dinkebach算法]

Source Code
Problem: 2728 User: orlando22
              Memory: 320K Time: 329MS//好菜的结果
Language: C++ Result: Accepted
  • Source Code
    #include <cmath>
    #include <iostream>
    
    //pku 2728 Desert king fulll graph->Prim
    using namespace std;
    
    const int MAXN=1001;
    
    //The length of the channel is the horizontal distance between the two villages.
    //The cost of the channel is the height of the lifter.
    struct Village{
        int x, y, z;
        double cost(const Village& _t){//the height of lifter
            return (int)fabs((double)(this->z-_t.z));
        }
        double len(const Village& _t){//the horizontal distance
            int dx=this->x-_t.x, dy=this->y-_t.y;
            return sqrt((double)(dx*dx+dy*dy));
        }
    };
    struct Edge{
        int i, j;
    };
    Village v[MAXN];
    bool inG[MAXN];
    int closest[MAXN], n;
    double lowcost[MAXN];
    Edge es[MAXN];
    
    //declaration
    void Prim(double rate);//complete graph is better than Kruskal
    double backtoEq();//using es[1...n] to calculate next rate
    
    //implementation
    void Prim(double rate){
        //array implementation without heap
        inG[1]=true;memset((void*)inG, 0, sizeof(bool)*(n+1));
        for(int i=2; i<n+1; i++){
            lowcost[i]=v[1].cost(v[i])-rate*v[1].len(v[i]), closest[i]=1, inG[i]=false;
        }
        //add for prim to get es[1...n-1]
        double minCost; int j;
        for(int i=1; i<n; i++){//add 1->n-1
            minCost=INT_MAX;
            //find the mimnal nodes
            for(int k=2; k<n+1; k++)
                if(lowcost[k]<minCost&&!inG[k])minCost=lowcost[k],j=k;
            inG[j]=true; es[i].i=j, es[i].j=closest[j];
            //update value
            for(int k=2; k<n+1; k++){
                double ncost=v[j].cost(v[k])-rate*v[j].len(v[k]);
                if(ncost<lowcost[k]&&!inG[k])lowcost[k]=ncost,closest[k]=j;
            }
        }
    }
    double backtoEq(){
        double dC=0, dL=0;
        for(int i=1; i<n; i++){
            dC+=v[es[i].i].cost(v[es[i].j]), dL+=v[es[i].i].len(v[es[i].j]);
        }
        return dC/dL;
    }
    
    int main()
    {
        //implementation
        while(cin>>n&&n!=0){
            for(int i=1; i<n+1; i++) cin>>v[i].x, cin>>v[i].y, cin>>v[i].z;
            //after initilation
            double nextRate=0, curRate;
            do{
                curRate=nextRate;
                Prim(curRate);
                nextRate=backtoEq();
            }while(fabs(curRate-nextRate)>1e-6);//sigh!
            printf("%.3lf\n",curRate);
        }
    
        return 0;
    }
    
Advertisements