Dinkelbach 算法程序源代码
谁有Dinkelbach 算法的程序源代码吗？ 对于0-1分数规划的Dinkelbach算法的分析武钢三中 吴豪[译]
摘要:
0-1分数规划问题是指求出解集{xi|xi=0或1}使目标(c1x1+c2x2+...+cnxn) /(d1x1+d2x2+…+dnxn)=cx/dx达到最大。对于分数规划问题，Dinkelbach提出了一个算法，它通过解决一个子问题Q(L)来得到原文题的解。这里Q是一个线性的最小化目标函数cx-Ldx，且满足x等于0或1。在本文中，我们证明了Dinkelbach算法在最坏情况下可以在O(log(nM))的时间内解决子问题，这里M=max{max|ci|,max|di|,1}。
1．0-1分数规划问题
要使两个线性函数的比值最大或最小的问题，我们称作分数规划问题或双曲线问题。分数规划问题在许多领域都可以找到。它在经济学中的应用有些常见的例子，如寻找最优收入比率或者在效益约束下的最佳物资调配问题。另外，系统效率也常常用比率来衡量，如收益/时间、利润/风险和消费/时间。有大量的文章对这类问题做了分析。
有几类分数规划问题已被广泛地研究。如0-1分数规划问题,它包含最优比率生成树问题，最优比率环问题，分数背包问题，以及分数剪枝问题。在本文中，我们研究0-1分数规划问题，它的描述如下：
令c=(c1,c2,…,cn)和d=(d1,d2,…,dn)为n维整数向量，那么一个0-1分数规划问题用公式描述如下:
FP: 最小化 (c1x1+…cnxn)/(d1x1…dnxn)=cx/dx
xi∈{0,1}
这里x表示列向量(x1,x2,…,xn)T .0-1值向量的子集Ω称作可行域，而x则是Ω的一个元素，我们称x为可行解。贯穿全文，我们假定对于任意可行解x,dx都是正数。这里我们记C=max{max|ci|,1},D=max{max|di|,1}。那么，显然问题的最优解在区间[-nC,nC]内。
对于分数规划问题，有许多算法都能利用下面的线性目标函数解决问题。
Q(L): 最小化 cx-Ldx
xi∈{0,1}
记z(L)为Q(L)的最值。令x*为分数规划的最优解，并且令L*= (cx*)/(dx*)(注：分数规划的最值)。那么下面就容易知道了：
z(L) > 0 当且仅当 L<L*
z(L) = 0 当且仅当 L=L*
z(L) < 0 当且仅当 L>L*
此外，Q(L*)的最优解也能使分数规划最优化。因此，解决分数规划问题在本质上等同于寻找L=L*使z(L)=0。出于这个目的，关于L的函数z(L)具有很多不错的性质：分段线性，凹函数，严格递减，z(-nC)<0，且z(nC)>0。根据上面的性质，显然当我们确定参量L，我们可以检验最值L*是否大于小于或等于当前的L。
有一些方法能够产生一系列收敛于L*的参量。其中一种借助于二分搜索。在两个不同的可行解的目标值不相同的情况下，他们的差距将大于等于1/(nD)^2。这暗示我们，当我们采用二分搜索时，最优值L*可以通过解决子问题Q(L)在最多O(log(2nC/(1/nD)^2))<=O(log(nCD))的时间内得到。
在1979年，Megiddo提出了一个巧妙的方法来系统地产生参量序列。他证明了如果子问题Q(L)能够通过O(p(n))的比较和O(q(n))的累加被解决，那么分数规划问题就能用O(p(n)+q(n))的时间被解决。
另一种方法理论上类似于牛顿迭代法，他被Isbell、Marlow和Dinkelbach提出（也被称作Dinkelbach算法）。这个算法在中被讨论（也可能是其他文献）。下一节将对它进行正式的论述。Schaible证明了对于非线性分数规划问题，二分搜索的方法的收敛速度仅仅是线性的，而Dinkelbach的收敛速度却是超线性的。另外，据说Dinkelbach算法在实际应用中强力而有效（参见的例子）。然而，Dinkelbach算法对于0-1分数规划问题的最坏时间复杂度却没有被证明。在本文中，我们证明了,Dinkelbach算法最多会在O(log(nCD))的时间内解决子问题。注意它的时间复杂度与普通的二分搜索相同。我们的结论暗示了，如果对于子问题Q(L)存在多项式算法，Dinkelbach算法也能够在多项式时间内解决分数规划问题。另外，即使子问题Q(L)是NP-完全或NP-难的，对于特殊的分数规划我们也能够在多项式时间内出解。
2．Dinkelbach算法的论述
它本质上是观察直线
z=cx’-Ldx’
于函数z(L)在L=L’处相切，这里x’是子问题Q(L’)的最优解。因此，-dx’是z(L)在L’处的斜率。而且很容易看出上面的直线与L轴相交与L=cx’/dx’.
现在我们来描述Dinkelbach对于分数规划的算法。Dinkelbach算法产生了收敛于L*的参量序列，如图1中细线所示的方式。
Dinkelbach算法：
步骤1：设L=L1,使 L*<=L1<=nC
步骤2：解决子问题Q(L)并得到最优解x
步骤3：如果z(L)=0，那么输出x并终止。否则，设L=cx/dx跳到步骤2
为了初始化L1,将用到nC，因此充分挖掘拓展问题的结构将能做出更好的选择。
3．Dinkelbach算法的分析
在这一节中，我们假定Dinkelbach算法终止于第k次。我们可以得到一个参量序列(L1,L2…Lk)和一个0-1值的向量(x1,x2…xk).z(L)的凸度暗示了下面的不等式：
dx1>dx2>…>dxk-1>=dxk>0
cx1+nCdx1>cx2+nCdx2>...>cxk-1+nCdxk-1>=cxk+nCxk>=0
由于函数z(L)是严格递减的，也很容易发现
z(L1)<z(L2)<...<z(Lk) = 0 并且 L1>L2>...>Lk
引理1 如果0>=z(Lr)>-1/nD (2<=r<=k) 那么 z(Lr)=0
证明 由于Lr=cxr-1/dxr-1
z(Lr)=cxr- Lr dxr=cxr-(cxr-1dxr)/(dxr-1)=(cxrdxr-1 – cxr-1dxr)/(dxr-1)
假定z(Lr)<0，那么(cxrdxr-1 – cxr-1dxr)<=-1。
因此不等式0<dxr-1<=nD暗示了z(Lr)<=-1/nD。它是矛盾的！
上面的引理来源于权向量c和d的完整性。这个引理暗示了如果z(Lr)>=-1/nD那么z(Lr)=0，因此该算法会中止于第r次。
引理2 如果0<=cxr+nCdxr<1,那么z(cxr/dxr)=0
证明 由于cxr+nCdxr是整数，如果0<=cxr+nCdxr<1，那么cxr+nCdxr=0并且cxr/dxr=-nC。由于最值L*>=-nC, xr是分数规划的最优解并且L*=cxr/dxr=-nC。那么显然有z(cxr/dxr)=z(L*)=0
上面的引理证明了如果cxr+nCdxr < 1，那么算法就在r或者r+1次终止。
现在给出主要引理：
引理3 如果1<=r<=k-1那么|z(Lr+1)|<=(1/2)|z(Lr)|或cxr+1+nCdxr+1 <=(1/2)(cxr+nCdxr)将满足。
证明 如果Lr+nC<=0 那么Lr=L*=-nC 并且它暗示着z(Lr)b=(1/2)z(Lr+1)=0
现在假定Lr+nC>0。就出现了两种情况：
情况(i) 首先我们来考虑z(Lr+1)( Lr+nC)<=z(Lr)/2*(Lr+1+nC)。这个情形如图２所示。在这个图中，直线z=cxr-Ldxr被记作lr。这里我们将用到图２的符号。令M为线段PR的中点。那么点M的坐标为( Lr+1 , z(Lr)( Lr+1+nC)/2/ Lr+nC )。因此条件z(Lr+1)( Lr+nC)<=z(Lr)*( Lr+1+nC)/2暗示着点Q=( Lr+1,z(Lr+1))在线段MR上。在这个条件下，我们证明不等式cxr+1+nC dxr+1<=(1/2)(cxr+nC, dxr)成立。这意味着直线lr+1与线段MR相交，lr+1也与线段M’R’相交，这里M’是线段P’R’的中点。现在我们证明这个不等式：
(Lr-Lr+1)( cxr+1+nC dxr+1)
=( cxr+1-Lr+1 dxr+1) (Lr+nC)- ( cxr+1-Lr dxr+1) (Lr+1+nC)
=z(Lr+1) (Lr+nC) - ( cxr+1-Lr dxr+1) (Lr+1+nC)
<= z(Lr) (Lr+1+nC)/2- ( cxr-Lr dxr) (Lr+1+nC)
=-(1/2) ( cxr-Lr dxr) (Lr+1+nC)=-(1/2)( cxr/ dxr - Lr) dxr( cxr/ dxr+nC)
=-(1/2)( Lr+1-Lr) ( cxr+nC dxr)=(1/2) (Lr-Lr+1) ( cxr+nC dxr)
由于Lr>Lr+1,那么不等式cxr+1+nC dxr+1<= (1/2) ( cxr+nC dxr)已经被证明。
情况(ii) 接着，考虑z(Lr+1)( Lr+nC)>z(Lr)/2*(Lr+1+nC)
|z(Lr+1)|=- z(Lr+1)< - z(Lr)(nC+ Lr+1)/2/(nC+ Lr)
=|z(Lr)|(1-(Lr- Lr+1)/ (nC+Lr))/2<=1/2 * |z(Lr)|
注意无论|z(L)|还是cx+nCdx在过程中都是不增长的。
在第一次，|z(L)|的值小于等于2n^2*CD。通过引理１，显然如果存在O(log(2n^2CD/(1/nD)))<=O(log(nCD))次，他们每个都能至少将z(L)减少50%那么，然后就能得到最优解。同样，引理２暗示了将cx+nCdx减少50%的次数O (log(2n^2CD))<=O(log(nCD))是最坏情况。引理３证明了每次将|z(L)|或cx+nCdx减少50%。因此，重复总次数的界限是O(log(nCD))。
定理４ Dinkelbanch算法最坏情况的运行次数是O(log(nCD))<=O(log(nM))，这里M=max{C,D}。
上面的定理证明了Dinkelbanch算法最坏运行次数是O(log(nCD))。它暗示了，如果对于Q(L)存在强多项式算法，Dinkelbach算法就能在多项式时间内解决分数规划问题。然而，当我们用多项式算法解决了子问题后，我们需要估计目标函数Q(L)的系数的输入规模。在下一节，我们将通过分析最优比率生成树和分数调配问题来讨论这一点。
4．讨论
Chandrasekaran提出了最优比率生成树的算法，它是基于二分搜索的。Dinkelbach算法可以在O(T(v,e)log(vCD))的时间解决该问题，这里v是点的个数，e是边的个数，并且用T(v,e)表示计算普通最小生成树的强多项式算法。很容易将Chandrasekaran的算法延伸到一般带有分数目标函数的矩阵胚规划问题。在这种情况下，在这种情形下，函数z(L)的断点数最大为n(n-1)/2（参见）因此，当可行域Ω是矩阵胚基础特征向量的集合。Dinkelbach算法就会在O(min{n^2,log(nCD)})后终止。
对于调配问题，已经研制了许多算法。大概最有名的算法就是Hungarian方法，并且它在最坏情况下的复杂度是O(v(v log v+e)) 。使用Hungarian方法，Dinkelbach算法可以用O(v(v log v+e)log(vCD))的时间解决分数调配问题。在中，Orlin和Ahuja提出了一个O(sqrt(v)e*log(vW))的算法来解决调配问题而且据说他们的算法因为强多项式算法而具有竞争性（参见也可）。在他们的算法中，它假定边权为整数，并且W代表边权的最大绝对值。为了将这个算法与Dinkelbach算法相结合，我们需要将在运行第r次解决的的子问题Q(L)用下面式子代替
(dxr-1) cx-( cxr-1) dx
这里xr-1表示代表运行第i次得到的最优解。因此，在每次运行中，Dinkelbach算法解决边权绝对值小于等于2v^2CD的调配问题。它暗示了Dinkelbach算法对于调配问题的最坏情况下的时间复杂度为
O(sqrt(v)e(log(2v^3CD))(log(eCD)))<=O(sqrt(v)e(log(vCD))^2)
我们说，如果一个具有线性目标函数的0-1整数规划问题存在多项式算法，我们可以利用上述目标函数在多项式时间内解决它的0-1分数规划问题。
致谢
作者在此感谢东京科学大学的Prof.Hirabayashi一直以来的鼓励与讨论。
0-1分数规划问题
分类： ACM•数学2012-05-01 20:01 307人阅读 评论(0) 收藏 举报
算法searchcini
0-1分数规划问题是指求出解集{xi|xi=0或1}使目标(c1x1+c2x2+...+cnxn) /(d1x1+d2x2+…+dnxn)=cx/dx达到最大。
对于分数规划问题，有许多算法都能利用下面的线性目标函数解决问题。
Q(L): 最小化 cx-Ldx xi∈{0,1}
记z(L)为Q(L)的最值。令x*为分数规划的最优解，并且令L*= (cx*)/(dx*)(注：分数规划的最值)。那么下面就容易知道了：
z(L) > 0 当且仅当 L<L*
z(L) = 0 当且仅当 L=L*
z(L) < 0 当且仅当 L>L*
因此，解决分数规划问题在本质上等同于寻找L=L*使z(L)=0。出于这个目的，关于L的函数z(L)具有很多不错的性质：分段线性，凹函数，严格递减，(-nC)<0，且z(nC)>0。
根据上面的性质，显然当我们确定参量L，我们可以检验最值L*是否大于小于或等于当前的L。二分搜索答案对于浮点误差控制较好，但是效率略低。
另一种方法类似于牛顿迭代的Dinkelbach算法，二分搜索的方法的收敛速度仅仅是线性的，而Dinkelbach的收敛速度却是超线性的。
Dinkelbach算法的论述
它本质上是观察直线
z=cx’-Ldx’
于函数z(L)在L=L’处相切，这里x’是子问题Q(L’)的最优解。因此，-dx’是z(L)在L’处的斜率。而且很容易看出上面的直线与L轴相交与L=cx’/dx’.
现在我们来描述Dinkelbach对于分数规划的算法。Dinkelbach算法产生了收敛于L*的参量序列，如图1中细线所示的方式。这里我们记C=max{max|ci|,1},显然问题的最优解在区间[-nC,nC]内。
Dinkelbach算法：
步骤1：设L=L1,使 L*<=L1<=nC
步骤2：解决子问题Q(L)并得到最优解x
步骤3：如果z(L)=0，那么输出x并终止。否则，设L=cx/dx跳到步骤2
以poj3111 K Best为例，选取k个元素，使得最大。
view plaincopy
double search(double p){
for(int i=0;i<n;i++)
jw.value=jw.v-jw.w*p;
Arrays.sort(jw);
double tv=0,tw=0;
for(int i=0;i<k;i++){
tv+=jw.v;
tw+=jw.w;
}
return tv/tw;
}
view plaincopy
double ans=1,temp;
while(true){
temp=search(ans);
if(Math.abs(ans-temp)<0.001)
break;
ans=temp;
}
0-1分数规划的应用：
1.最优比率生成树（poj2526 Desert King）：一个带权完全图，每条边都有自己的花费值cost和 长度dis，求一个生成树，使得r=（∑cost*x ） / （∑dis*x ）最小。
解法：search时将每条边边权设为cost-dis*p，求最小生成树时选取的边既为一组解，不断迭代。
2.最优比率环：(poj3621 SightSeeing)
有向图中，每个点有一个权值C(1<=c<=1000)，每条边有权值D(1<=D<=1000 要求找出一个环（没有重点）， 环包含的点数要>=2，使得sum{C}/sum{D}最大。
解法：由于此题求解的是存在性问题，不可套用Dinkelbach算法迭代，因此二分答案求解最大值，对于给定的g如果以<a,b>权重为g*D-C的图中存在环路，则ans>=g（注意二分double时结束条件right-left<eps,转移为left=mid、right=mid）。
3.最大密度子图
详见《最大权闭合图&&最大密度子图》
01分数规划的理解以及练习poj2976 poj 3621 poj 2728
分类： 数学 poj2012-10-08 14:21 282人阅读 评论(0) 收藏 举报
distancestructc算法图形编程
01分数规划的思想的描述如下：令c=(c1,c2,…,cn)和d=(d1,d2,…,dn)为n维整数向量，那么一个0-1分数规划问题用公式描述如下:FP: 最小化或者最大化(c1x1+…cnxn)/(d1x1…dnxn)=cx/dx xi∈{0,1}。这个问题的解决已经有了成形的算法，但是网上的解释大多过难难以直接理解，大部分对于数学基础的要求是很高的。
先说说这种思想能够解决哪些问题 1：裸的01分数规划pku2976适合刚刚学习的练手以及试模板 2：01分数规划与图形的结合比如说环的判断pku3621比如最优比例生成树pku2728 解决的思想大致一样 关键是对结论的简要证明，加深理解，然后灵活的运用。
下面谈谈解决01分数规划的方法以及结论
首先，假设要c*x/d*x最大（x，c，d均指向量），假设最优解为 val，设置这样的子问题，z(L)=c*x-d*x*L 分配x向量使 z最大，那么可以知道 如果c*x-L*x*d<0 那么无论怎么分配x向量都无法满足L值即val必然小于L 同样 我们看大于0的情况 这个时候说明可以分配x使得目标解大于L那么val必然大于L，那么很显然 L=val的时候c*x-d*x*L=0，并且z的值具有单调性下面是证明
设L1,L2,且L1>L2
z(L1) = max { c*x-L1*d*x } = c*x1-L1*d*x1(x1是使得其最大的向量) <c*x1-L2*d*x1≤ max{ c*x-d*x*L2 } = z(L2)
=> z(l)是关于L的单调递减函数。
因此我们可以用二分或者迭代来解决
这里只是简单的证明并没有很严格的证明 但是对于使用来讲已经足够。
下面具体使用结论来解决上面的三道例子；
pku2976 很直接的一道题 题意就是简单01分数规划
来说说做法 由于题目已经告诉我们去掉K个x，我们要使得分数最高，设置子问题 z,根据上面的公式很容易写出程序。
view plaincopy
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cstdio>
#include<vector>
#include<sstream>
#include<string>
#include<climits>
#include<stack>
#include<set>
#include<bitset>
#include<cmath>
#include<deque>
#include<map>
#include<queue>
#define iinf 2000000000
#define linf 1000000000000000000LL
#define dinf 1e200
#define eps 1e-6
#define all(v) (v).begin(),(v).end()
#define sz(x) x.size()
#define pb push_back
#define mp make_pair
#define lng long long
#define sqr(a) ((a)*(a))
#define pii pair<int,int>
#define pll pair<lng,lng>
#define pss pair<string,string>
#define pdd pair<double,double>
#define X first
#define Y second
#define pi 3.14159265359
#define ff(i,xi,n) for(int i=xi;i<=(int)(n);++i)
#define ffd(i,xi,n) for(int i=xi;i>=(int)(n);--i)
#define ffl(i,r) for(int i=head;i!=-1;i=edge.next)
#define cc(i,j) memset(i,j,sizeof(i))
#define two(x) ((lng)1<<(x))
#define N 555555
#define M 1000000
#define lson l , mid , rt << 1
#define rson mid + 1 , r , rt << 1 | 1
#define Mod n
#define Pmod(x) (x%Mod+Mod)%Mod
using namespace std;
typedef vector<int> vi;
typedef vector<string> vs;
typedef unsigned int uint;
typedef unsigned lng ulng;
template<class T> inline void checkmax(T &x,T y)
{
if(x<y) x=y;
}
template<class T> inline void checkmin(T &x,T y)
{
if(x>y) x=y;
}
template<class T> inline T Min(T x,T y)
{
return (x>y?y:x);
}
template<class T> inline T Max(T x,T y)
{
return (x<y?y:x);
}
template<class T> T gcd(T a,T b)
{
return (a%b)==0?b:gcd(b,a%b);
}
template<class T> T lcm(T a,T b)
{
return a*b/gcd(a,b);
}
template<class T> T Abs(T a)
{
return a>0?a:(-a);
}
using namespace std;
struct node
{
int x,y;
} a;
double low,high,mid;
int cmp(node bb,node cc)
{
return mid*bb.y-bb.x<mid*cc.y-cc.x;
}
int main()
{
int n,k;
double sub,par,res,ans;
while(scanf("%d%d",&n,&k)!=EOF)
{
if(n==0&&k==0) break;
int i,j;
for(i=0; i<n; i++) scanf("%d",&a.x);
for(i=0; i<n; i++) scanf("%d",&a.y);
low=0;
high=1000;
ans=0;
while(low+eps<high)
{
mid=(low+high)/2.0;
sort(a,a+n,cmp);
sub=par=0;
for(i=0; i<n-k; i++)
{
sub+=a.x;
par+=a.y;
}
ans=sub/par;
res=mid*par-sub;
if(res<0) low=mid;
else high=mid;
}
printf("%.0lf\n",100*ans);
}
return 0;
}
pku3621 这道题是要求一个最大（边权和/点权和)的环 同样设置子问题 边权-点权*L最大的回路，如果存在并大于0 r=mid 可以直接用spfa的判环 将权编程负值然后判断有没有负环即可
view plaincopy
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cstdio>
#include<vector>
#include<sstream>
#include<string>
#include<climits>
#include<stack>
#include<set>
#include<bitset>
#include<cmath>
#include<deque>
#include<map>
#include<queue>
#define iinf 2000000000
#define linf 1000000000000000000LL
#define dinf 1e200
#define eps 1e-3
#define all(v) (v).begin(),(v).end()
#define sz(x) x.size()
#define pb push_back
#define mp make_pair
#define lng long long
#define sqr(a) ((a)*(a))
#define pii pair<int,int>
#define pll pair<lng,lng>
#define pss pair<string,string>
#define pdd pair<double,double>
#define X first
#define Y second
#define pi 3.14159265359
#define ff(i,xi,n) for(int i=xi;i<=(int)(n);++i)
#define ffd(i,xi,n) for(int i=xi;i>=(int)(n);--i)
#define ffl(i,r) for(int i=head;i!=-1;i=edge.next)
#define cc(i,j) memset(i,j,sizeof(i))
#define two(x) ((lng)1<<(x))
#define N 1050
#define M 10000
#define lson l , mid , rt << 1
#define rson mid + 1 , r , rt << 1 | 1
#define Mod n
#define Pmod(x) (x%Mod+Mod)%Mod
using namespace std;
typedef vector<int> vi;
typedef vector<string> vs;
typedef unsigned int uint;
typedef unsigned lng ulng;
template<class T> inline void checkmax(T &x,T y){if(x<y) x=y;}
template<class T> inline void checkmin(T &x,T y){if(x>y) x=y;}
template<class T> inline T Min(T x,T y){return (x>y?y:x);}
template<class T> inline T Max(T x,T y){return (x<y?y:x);}
template<class T> T gcd(T a,T b){return (a%b)==0?b:gcd(b,a%b);}
template<class T> T lcm(T a,T b){return a*b/gcd(a,b);}
template<class T> T Abs(T a){return a>0?a:(-a);}
struct pp
{
int v,next,w;
}edge;
double dist;
int f,head;
bool inque;
int cnt;
int n,m,tot;
double l,r,ans,mid;
inline void addedge(int u,int v,int w)
{
edge.v=v,edge.w=w,edge.next=head,head=tot++;
}
bool spfa()
{
queue<int> q;
ff(i,1,n) dist=iinf;
cc(cnt,0);
cc(inque,0);
q.push(1);
inque=1,cnt=1;
dist=0;
while(!q.empty())
{
int u=q.front();
inque=0;
q.pop();
ffl(i,u)
{
double w=edge.w;
int v=edge.v;
w=mid*w-f;
if(dist>dist+w)
{
dist=dist+w;
if(!inque)
{
inque=1;
q.push(v);
cnt++;
if(cnt>=n) return 1;
}
}
}
}
return 0;
}
int main()
{
while(scanf("%d%d",&n,&m)==2)
{
ff(i,1,n) scanf("%d",f+i);
tot=0;
cc(head,-1);
ff(i,1,m)
{
int u,v,w;
scanf("%d%d%d",&u,&v,&w);
addedge(u,v,w);
}
l=0,r=1996;
while(r-l>=0.001)
{
mid=(l+r)/2;
if(spfa())
{
l=mid;
ans=mid;
}
else
r=mid;
}
printf("%.2f\n",ans);
}
return 0;
}
pku2728
最优比例生成树
说白了就是求最小的(代价权和/边权和)的生成树问题 同样设置子问题 代价权-边权*L最小 即重新构造图 然后求最小生成树 同样如果 小于0 自己推导下就会知道 r=mid 另外poj会卡stl 裸奔的prim算法能过 优先队列优化的反而超时
view plaincopy
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cstdio>
#include<vector>
#include<sstream>
#include<string>
#include<climits>
#include<stack>
#include<set>
#include<bitset>
#include<cmath>
#include<deque>
#include<map>
#include<queue>
#define iinf 2000000000
#define linf 1000000000000000000LL
#define dinf 1e200
#define eps 1e-5
#define all(v) (v).begin(),(v).end()
#define sz(x) x.size()
#define pb push_back
#define mp make_pair
#define lng long long
#define sqr(a) ((a)*(a))
#define pii pair<int,int>
#define pll pair<lng,lng>
#define pss pair<string,string>
#define pdd pair<double,double>
#define X first
#define Y second
#define pi 3.14159265359
#define ff(i,xi,n) for(int i=xi;i<=(int)(n);++i)
#define ffd(i,xi,n) for(int i=xi;i>=(int)(n);--i)
#define ffl(i,r) for(int i=head;i!=-1;i=edge.next)
#define cc(i,j) memset(i,j,sizeof(i))
#define two(x) ((lng)1<<(x))
#define N 1050
#define M 1000000
#define lson l , mid , rt << 1
#define rson mid + 1 , r , rt << 1 | 1
#define Mod n
#define Pmod(x) (x%Mod+Mod)%Mod
using namespace std;
typedef vector<int> vi;
typedef vector<string> vs;
typedef unsigned int uint;
typedef unsigned lng ulng;
template<class T> inline void checkmax(T &x,T y){if(x<y) x=y;}
template<class T> inline void checkmin(T &x,T y){if(x>y) x=y;}
template<class T> inline T Min(T x,T y){return (x>y?y:x);}
template<class T> inline T Max(T x,T y){return (x<y?y:x);}
template<class T> T gcd(T a,T b){return (a%b)==0?b:gcd(b,a%b);}
template<class T> T lcm(T a,T b){return a*b/gcd(a,b);}
template<class T> T Abs(T a){return a>0?a:(-a);}
int n,h;
struct ps
{
int x,y;
}cor;
inline double Distance(ps &p, ps &q )
{
return sqrt((double)sqr(p.x-q.x)+(double)sqr(p.y-q.y));
}
double l,r,ans,mid,dist;
double cost,len,ww;
inline void prim()
{
ans=0;
bool in={};
dist=0;
ff(i,2,n) dist=iinf;
ff(i,1,n)
{
double w=iinf;
int id;
ff(j,1,n)
if(!in&&dist<w)
w=dist,id=j;
in=1;
ans+=w;
ff(j,1,n)
if(!in&&ww<dist) dist=ww;
}
}
int main()
{
while(scanf("%d",&n)==1&&n)
{
ff(i,1,n)
scanf("%d%d%d",&cor.x,&cor.y,&h);
ff(i,1,n)
ff(j,i+1,n)
cost=cost=Abs(h-h),len=len=Distance(cor,cor);
l=0,r=500000;
while(r-l>=eps)
{
mid=(l+r)/2;
for(int i=1;i<=n;i++)
for(int j=i+1;j<=n;j++)
ww=ww=cost-len*mid;
prim();
if(ans<0) r=mid;
else
l=mid;
}
printf("%.3f\n",mid);
}
return 0;
}
这类问题的种类很少 基本属于模板题范畴 不如动态规划这些博大精深···
Poj 2728 最优比率生成树
/*
目标:min{∑costi/∑leni}
逼近的思想，∑costi/∑leni<=x，即 ∑(costi-x*leni)<=0 是一个单调递减函数
即求边为costi-x*leni的 MST
可以用二分，但比较慢
用迭代快好多
*/
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<cstdlib>
using namespace std;
const double esp=0.00001;
const int MAXN=1010;
const double DINF=1000000000.0;
struct Point{
int x,y,z;
}points;
int N;
bool vi;
double dist;
int pre;
double cal(int a,int b){
return sqrt(1.*(points.x-points.x)*(points.x-points.x)+
1.*(points.y-points.y)*(points.y-points.y));
}
double prim(double x){
memset(vi,0,sizeof(vi));
for(int i=2;i<=N;i++){
dist=abs(points.z-points.z)-cal(1,i)*x;
pre=1;
}
dist=0;vi=1;
double cost=0,len=0;
for(int i=1;i<N;i++){
double Min=DINF;
int u;
for(int j=2;j<=N;j++)
if(!vi&&Min>dist){
Min=dist;
u=j;
}
vi=1;
cost+=abs(points].z-points.z);
len+=cal(pre,u);
for(int j=2;j<=N;j++){
double val=abs(points.z-points.z)-cal(u,j)*x;
if(!vi&&dist>val){
dist=val;
pre=u;
}
}
}
return cost/len;
}
int main(){
while(scanf("%d",&N),N){
for(int i=1;i<=N;i++)
scanf("%d%d%d",&points.x,&points.y,&points.z);
//分数规划用：Dinkelbach算法
//每次迭代子问题的解cost`/len`进去，这样会不断逼近最优解
double a=0,b;
while(1){
b=prim(a);
if(fabs(b-a)<esp)break;
a=b;
}
printf("%.3f\n",b);
}
return 0;
}
/*
//cost-len*x<=0
double low=0,high=100.0; //其实二分20多次已经很足够了
while(high-low>esp){
double mid=(low+high)/2;
if(prim(mid))high=mid;
else low=mid;
}
printf("%.3f\n",high);
*/
页:
[1]