【Learning】分数规划

分数规划

　　
​　　分数规划是一类决策性问题。一般地，题目会要求你针对问题规划一种方案，使得其代价函数最小或最大。其中，代价函数一般是分数形式，且分子分母的构成元素一般呈现一一对应关系。
　　
　　直接上例题观察：BZOJ2402
　　
　　
　　
　　分数规划的重要思路是二分答案。我们记答案为(c)，则问题转变为判断是否存在((i,j))使得

[frac{y_i+q_j}{x_i+p_j}ge c\ y_i+q_jge c(x_i+p_j)\ (y_i-c*x_i)+(q_j-c*p_j)ge0 ]

​　　化简后我们发现，原本分别处于分子和分母的元素有了一一对应的关系，使得每一类元素之间在表达式上互不干扰。只要知道此式有没有可能成立，就能判断(c)是否合法。这也是分数规划的一个重要突破点。
　　
​　　所以显然是要分别找到(a)到(b)的路径上两个点(i,j)，使得((y_i-c*x_i))和((q_j-c*p_j))分别取最大值。如果两者加起来大于等于0，说明(c)这个答案合法，因此我们提升二分下界；否则(c)答案不合法，我们降低二分上界。
　　
　　取最大值的这些操作，可以用树剖+上凸包回答。
　　
　　总的时间复杂度是(O(nlog^4n))，但是能过！
　　

Code

　　

#include <cstdio>
#include <vector>
#include <algorithm>
#define mp make_pair
#define pb push_back
using namespace std;
typedef pair<double,double> pdd;
const int N=30005;
const double EPS=1e-6;
int n;
double a[N][4];
int h[N],tot,dep[N],pre[N][16],size[N],son[N];
struct Edge{int v,next;}e[N*2];
int who[N],loc[N],locnt,top[N];
inline void addEdge(int u,int v){
	e[++tot]=(Edge){v,h[u]}; h[u]=tot;
	e[++tot]=(Edge){u,h[v]}; h[v]=tot;
}
void dfs1(int u,int fa){
	size[u]=1; son[u]=0;
	dep[u]=dep[fa]+1;
	pre[u][0]=fa;
	for(int i=1;i<=15;i++) pre[u][i]=pre[pre[u][i-1]][i-1];
	for(int i=h[u],v;i;i=e[i].next)
		if((v=e[i].v)!=fa){
			dfs1(v,u);
			size[u]+=size[v];
			if(!son[u]||size[v]>size[son[u]]) son[u]=v;
		}
}
void dfs2(int u,int _top){
	top[u]=_top;
	who[loc[u]=++locnt]=u;
	if(!son[u]) return;
	dfs2(son[u],_top);
	for(int i=h[u],v;i;i=e[i].next)
		if((v=e[i].v)!=pre[u][0]&&v!=son[u])
			dfs2(v,v);
}
inline void merge(pdd &u,pdd v){
	if(v.first>u.first) u.first=v.first;
	if(v.second>u.second) u.second=v.second;
}
inline double getk(pdd u,pdd v){return (v.second-u.second)/(v.first-u.first);}
struct Hull{
	vector<pdd> a;
	void insert(double x,double y){a.pb(mp(x,y));}
	void build(){
		sort(a.begin(),a.end());
		int top=0;
		vector<pdd> st;
		for(int i=0,sz=a.size();i<sz;i++){
			while(top>1&&getk(st[top-2],st[top-1])<getk(st[top-1],a[i]))
				top--,st.pop_back();
			st.pb(a[i]); top++;
		}
		a=st;
		a.resize(top);
	}
	double query(double k){
		int l=0,r=a.size()-2,mid;
		while(l<=r){
			mid=(l+r)>>1;	
			if(getk(a[mid],a[mid+1])-EPS<=k) 
				r=mid-1;
			else l=mid+1;
		}
		return a[l].second-k*a[l].first;
	}
};
namespace SEG{
	int rt,sz,ch[N*2][2];
	Hull s[N*2][2];
	void build(int &u,int l,int r){
		u=++sz;
		for(int i=l;i<=r;i++){
			int x=who[i];
			s[u][0].insert(a[x][0],a[x][1]);
			s[u][1].insert(a[x][2],a[x][3]);
		}
		s[u][0].build(); s[u][1].build();
		if(l==r) return;
		int mid=(l+r)>>1;
		build(ch[u][0],l,mid);
		build(ch[u][1],mid+1,r);
	}
	pdd query(int u,int l,int r,int L,int R,double k){
		if(L<=l&&r<=R)
			return mp(s[u][0].query(k),s[u][1].query(k));
		pdd res=mp(-1e10,-1e10);	
		int mid=(l+r)>>1;
		if(L<=mid) merge(res,query(ch[u][0],l,mid,L,R,k));
		if(mid<R) merge(res,query(ch[u][1],mid+1,r,L,R,k));
		return res;
	}
}
int getLCA(int u,int v){
	if(dep[u]<dep[v]) swap(u,v);
	for(int i=15;i>=0;i--)
		if(dep[pre[u][i]]>=dep[v]) u=pre[u][i];
	if(u==v) return u;
	for(int i=15;i>=0;i--)
		if(pre[u][i]!=pre[v][i]) u=pre[u][i],v=pre[v][i];
	return pre[u][0];
}
bool judge(int u,int v,double k){
	pdd res=mp(-1e10,-1e10);	
	int lca=getLCA(u,v);		
	for(;dep[top[u]]>=dep[lca];u=pre[top[u]][0])
		merge(res,SEG::query(SEG::rt,1,n,loc[top[u]],loc[u],k));
	if(dep[u]>=dep[lca])
		merge(res,SEG::query(SEG::rt,1,n,loc[lca],loc[u],k));
	for(;dep[top[v]]>=dep[lca];v=pre[top[v]][0])
		merge(res,SEG::query(SEG::rt,1,n,loc[top[v]],loc[v],k));
	if(dep[v]>=dep[lca])
		merge(res,SEG::query(SEG::rt,1,n,loc[lca],loc[v],k));
	return res.first+res.second+EPS>=0;	
}
int main(){
	scanf("%d",&n);
	for(int i=0;i<4;i++)
		for(int j=1;j<=n;j++) scanf("%lf",&a[j][i]);
	for(int i=1,u,v;i<n;i++){
		scanf("%d%d",&u,&v);
		addEdge(u,v);
	}
	dfs1(1,0);
	dfs2(1,1);
	SEG::build(SEG::rt,1,n);
	int m,u,v;
	scanf("%d",&m);
	while(m--){
		scanf("%d%d",&u,&v);
		double l=0,r=100001,mid,eps=5e-4;
		while(l+eps<r){
			mid=(l+r)*0.5;
			if(judge(u,v,mid)) l=mid;
			else r=mid;
		}
		printf("%.5lf
",l);
	}
	return 0;
}

　　
　　　　
　　

例题

　　

【BZOJ4819】【SDOI2017】新生舞会

　　
​　　传送门
　　
　　按照分数规划的基本思想，我们二分答案(c)，化简题目中的代价函数式，转化成一个判定问题：

[frac{x_1+x_2+...+x_n}{y_1+y_2+...+y_n}ge c\ (x_1-c*y_1)+(x_2-c*y_2)+...+(x_n-c*y_n)ge 0 ]

　　其中(x_i)表示第(i)个男生和他所选女生的喜悦程度，(y_i)表示对应的不协调程度。
　　
​　　我们发现虽然各个括号在表达式上互不相关，但是从题意的角度来看：不能有多个男生同时选定一个女生。
　　
​　　但是我们的目的还是要将这(n)个括号的和最大化，并和0比较，从而判定答案(c)是否合法。
　　
​　　匹配，带权，最大化——我们想到了最大带权匹配模型。可以从最大费用最大流的角度来建图：源点向所有男生连((1,0))的边；所有男生向所有女生连边，边权为相应的括号，流量为1；所有女生向汇点连((1,0))的边。
　　
​　　那么我们跑一次费用流就可以得出一种方案，使得括号之和最大。此时和0判断即可。当然用KM会快很多。
　　　

Code

　　　

#include <cstdio>
#include <queue>
using namespace std;
const int N=105,INF=1e9;
const double EPS=1e-6;
int n,a[N][N],b[N][N];
int h[N*2],tot;
struct Edge{int v,f;double c;int next;}e[(N*N+2*N)*2];
int S,T;
double dis[N*2];
int pree[N*2],preu[N*2];
bool inq[N*2];
queue<int> q;
inline void addEdge(int u,int v,int f,double c){
    e[++tot]=(Edge){v,f,c,h[u]}; h[u]=tot;
    e[++tot]=(Edge){u,0,-c,h[v]}; h[v]=tot;
}
bool spfa(){
    for(int i=1;i<=T;i++)
        dis[i]=-INF,inq[i]=false;
    while(!q.empty()) q.pop();
    q.push(S);
    dis[S]=0; inq[S]=true;
    while(!q.empty()){
        int u=q.front(); q.pop();
        inq[u]=false;
        for(int i=h[u],v;i;i=e[i].next)
            if(e[i].f&&dis[v=e[i].v]<dis[u]+e[i].c){
                dis[v]=dis[u]+e[i].c;
                preu[v]=u; pree[v]=i;
                if(!inq[v]){
                    inq[v]=true;
                    q.push(v);
                }
            }
    }
    return dis[T]>-INF;
}
double maxcostflow(){
    double res=0;
    while(spfa()){
        int flow=INF;
        for(int u=T;u!=S;u=preu[u])
            flow=min(flow,e[pree[u]].f);
        res+=dis[T]*flow;
        for(int u=T;u!=S;u=preu[u]){
            e[pree[u]].f-=flow;
            e[pree[u]^1].f+=flow;
        }
    }
    return res;
}
bool judge(double c){
    S=n*2+1; T=S+1;
    for(int i=1;i<=T;i++) h[i]=0;
    tot=1;
    for(int i=1;i<=n;i++)
        addEdge(S,i,1,0),addEdge(n+i,T,1,0);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            addEdge(i,n+j,1,a[i][j]-c*b[i][j]);
    return maxcostflow()+EPS>=0;
}
int main(){
    scanf("%d",&n);
    for(int i=1;i<=n;i++) 
        for(int j=1;j<=n;j++) scanf("%d",&a[i][j]);
    for(int i=1;i<=n;i++) 
        for(int j=1;j<=n;j++) scanf("%d",&b[i][j]);
    double l=0,r=10000,mid;
    while(l+(5e-8)<r){
        mid=(l+r)*0.5;
        if(judge(mid)) l=mid;
        else r=mid;
    }
    printf("%.6lf
",l);
    return 0;
}

　　
​　　
　　
​　　
　　

【BZOJ1690】【Usaco2007 Dec】奶牛的旅行

　　
　　传送门（权限题）

　　按照套路二分答案展开化简，得到一个用于判定(c)是否合法的式子：

[(a_1-c*b_1)+(a_2-c*b_2)+...+(a_{m-1}-c*b_{m-1})ge 0 ]

　　其中选择的路径有(m)个点，(a_i)表示第(i)个点的乐趣值，(b_i)表示第(i)个点走到第(i+1)个点的距离。由于第(m)个点和第1个点相同，价值又不重复计算，所以(a_m)不需要算进答案，式子变得十分规律好看。
　　
　　那么我们把每一条边((u,v,w))的边权设置为(a_u-c*w)。如果图中存在正环，那么选择这一个环游览就能满足上述式子。
　　
　　所以用spfa来判断每次二分判定时构建的图是否存在正环即可。
　　　

Code

　　

#include <cstdio>
#include <queue>
using namespace std;
const int N=1005,M=5005,INF=1e9;
const double EPS=1e-6;
int n,m,a[N],inp[M][3];
int h[N],tot;
struct Edge{int v;double w;int next;}e[M];
inline void addEdge(int u,int v,double w){e[++tot]=(Edge){v,w,h[u]}; h[u]=tot;}
double dis[N];
bool inq[N];
int qt[N];
queue<int> q;
bool spfa(){
    while(!q.empty()) q.pop();
    for(int i=1;i<=n;i++) dis[i]=0,inq[i]=true,qt[i]=1,q.push(i);
    while(!q.empty()){
        int u=q.front(); q.pop();
        inq[u]=false;
        for(int i=h[u],v;i;i=e[i].next)
            if(dis[v=e[i].v]+EPS>dis[u]+e[i].w){
                dis[v]=dis[u]+e[i].w;
                if(!inq[v]){
                    inq[v]=true;
                    qt[v]++;
                    if(qt[v]>n) return true;
                    q.push(v);
                }
            }
    }
    return false;
}
bool judge(double c){
    tot=0;
    for(int i=1;i<=n;i++) h[i]=0;
    for(int i=1;i<=m;i++)
        addEdge(inp[i][0],inp[i][1],-(a[inp[i][0]]-c*inp[i][2]));
    return spfa();
}
int main(){
    scanf("%d%d",&n,&m);
    for(int i=1;i<=n;i++) scanf("%d",&a[i]);
    for(int i=1;i<=m;i++) scanf("%d%d%d",inp[i],inp[i]+1,inp[i]+2);
    double l=0,r=1000,mid;
    while(l+(1e-3)<r){
        mid=(l+r)*0.5;
        if(judge(mid)) l=mid;
        else r=mid;
    }
    printf("%.2lf
",l);
    return 0;
}

　　
　　
　　
　　
　　

【BZOJ4898】【Apio2017】商旅

　　
​　　传送门
　　
​　　此题和上题非常类似。但是题目的“环”好像没有强调是简单环，只需要走一圈回到出发点即可。
　　
​　　首先用Floyd计算出原图上两两点之间的距离(w_{u,v})。
　　
​　　再考虑购买商品的过程。题目所谓的“从A地买一个商品，拿着走向B地卖掉，利润为x”，可以简化为从A向B连一条权值为x的边。显然从A走向B总是应该选利润最大的那个商品走，所以A向B只需要连一条边。
　　
​　　当然，为了满足题目意思，对于原图的每一条边我们应该相应地连一条两端一样，边权为0的边，以适应空手旅行的情况。
　　
​　　原问题变成：在新构建的图上，找到一个环，使得代价函数最大。
　　
　　惯例化简代价函数，得到判定式

[(x_1-c*w_1)+(x_2-c*w_2)+...+(x_{m-1}-c*w_{m-1})ge 0 ]

　　(x_i)表示第(i)个点走向第(i+1)个点的边权（也就是这么走的利润），(w_i)表示从(i)走到(i+1)的距离（和上一题几乎一模一样）。
　　
​　　所以再次用spfa判正环是否存在即可。
　　
　　注意题目要求输出整数，我们用整数进行二分，可以提高效率。
　　

Code

　　

#include <cstdio>
using namespace std;
const int N=105,M=9905,K=1005,INF=1e9;
const double EPS=1e-6;
int n,m,kk,a[N],w[N][N];
int sell[N][K],buy[N][K];
int cnt,inp[N*N*2][3];
int h[N],tot;
struct Edge{int v;double w;int next;}e[N*N*2];
inline void addEdge(int u,int v,double w){e[++tot]=(Edge){v,w,h[u]}; h[u]=tot;}
double dis[N];
bool inq[N];
inline int min(int x,int y){return x<y?x:y;}
inline int max(int x,int y){return x>y?x:y;}
bool dfs(int u){
    inq[u]=true;
    for(int i=h[u],v;i;i=e[i].next){
        v=e[i].v;
        if(dis[u]+e[i].w<dis[v]+EPS){
            dis[v]=dis[u]+e[i].w;
            if(inq[v]) return true;
            else if(dfs(v)) return true;
        }
    }
    inq[u]=false;
    return false;
}
bool judge(double c){
    tot=0;
    for(int i=1;i<=n;i++) h[i]=0;
    for(int i=1;i<=cnt;i++)
        addEdge(inp[i][0],inp[i][1],-(inp[i][2]-c*w[inp[i][0]][inp[i][1]]));
    for(int i=1;i<=n;i++) dis[i]=0,inq[i]=false;
    for(int i=1;i<=n;i++) 
        if(dfs(i)) return true;
    return false;
}
int main(){
    scanf("%d%d%d",&n,&m,&kk);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=kk;j++) scanf("%d%d",&buy[i][j],&sell[i][j]);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++) w[i][j]=INF;
    for(int i=1,u,v,c;i<=m;i++){
        scanf("%d%d%d",&u,&v,&c);
        w[u][v]=min(w[u][v],c);
        cnt++;
        inp[cnt][0]=u; inp[cnt][1]=v; inp[cnt][2]=0;
    }
    for(int k=1;k<=n;k++)
        for(int i=1;i<=n;i++)
            for(int j=1;j<=n;j++)
                w[i][j]=min(w[i][j],w[i][k]+w[k][j]);
    for(int u=1;u<=n;u++)
        for(int v=1;v<=n;v++)
            if(w[u][v]!=INF){
                int best=0;
                for(int i=1;i<=kk;i++)
                    if(buy[u][i]!=-1&&sell[v][i]!=-1)
                        best=max(best,sell[v][i]-buy[u][i]);
                if(best){
                    cnt++;
                    inp[cnt][0]=u; inp[cnt][1]=v; inp[cnt][2]=best;
                }
            }
    int l=0,r=1e9,mid;
    while(l<=r){
        mid=(l+r)>>1;
        if(judge(mid)) l=mid+1;
        else r=mid-1;
    }
    printf("%d
",r);
    return 0;
}

　　
　　
　　
　　　　
　　

细节注意

　　
　　题目要求的精度非常关键，不要因为二分的控制部分不精确，使得答案出现偏差。
　　
　　一定要记得使用(eps)！比如在spfa判定是否要松弛的地方要加上eps。