【BZOJ 1701】Cow School（斜率优化/动态凸包/分治优化）

原题题解和数据下载 Usaco2007 Jan

题意

小牛参加了n个测试，第i个测试满分是(p_i),它的得分是(t_i)。老师去掉(t_i/p_i)最小的d个测试，将剩下的总得分/总满分作为小牛的得分。

小牛想知道多少个d存在比老师计算的分数更高的选择测试的方案，并输出这些d。

题解

基础思路

排好序后，$ frac {t_1} {p_1} < frac {t_2} {p_2}<..< frac {t_n} {p_n}$。

如果d==j，老师给的分数是(r_j=frac {t_{j+1}+t_{j+2}+..t_{n}} {p_{j+1}+p_{j+2}+..p_{n}}=frac{st_j}{sp_j})。

要找的就是是否存在一个集合(S(|S|=n-j))，满足：

[frac {sum_{iin S}t_i} {sum_{iin S}p_i}>r_j ]

等价于(sum_{iin S}t_i- {sum_{iin S}p_i} cdot r_j>0)。也就是：

[sum_{iin S}(t_i-p_i cdot r_j)>0 Rightarrow sum_{iin S}t_icdot sp_j-p_icdot st_j>0 ]

贪心地找最大的 n-j 个加起来，复杂度是(O(n^2log n))。

法1 平衡树维护动态凸包

一个优化的方法是，考虑(z=t_icdot sp_j-p_icdot st_j)在(1..j)里最大值(g[j]) 和(j+1..n)里的最小值(f[j])。

如果(f[j]) 小于(g[j])，那么意味着从(j+1..n)中取出最小值换为(1..j)中的最大值，可以更优。

(g[j]=max{a_jcdot x_i+b_jcdot y_i})实际上就是斜率优化的形式。

这里(x_i=p_i,y_i=t_i,a_j=-st_j,b_j=sp_j)。

直线为(y=st_j/sp_jcdot x+z/sp_j)。

所以每个点坐标为((p_i,t_i))，要找一个点，斜率为(-a_j/b_j)的直线经过它时，纵截距最大。需要维护一个上凸壳。x不是单调的，需要用平衡树来维护这个凸壳。

(f[j])同理求。总复杂度(O(nlog n))。

法2 普通维护凸包

实际上还可以更优。

求(g[j])时加入上凸壳的点((p_j,t_j))和原点连线的斜率(t_j/p_j)，一定比之前加入的任意点((p_{j'},t_{j'}))的要大，于是新点一定可以保留在凸壳上。直线斜率(st_j/sp_j)一定比(t_j/p_j)大。于是当前点右边的点就没有用了。因为直线斜率递增，所以最优解的位置递减。

求(f[j])时加入下凸壳的点((p_j,t_j))和原点连线的斜率(t_j/p_j)，一定比之前加入的任意点((p_{j'},t_{j'}))的要小，于是新点一定可以保留在凸壳上。直线斜率(st_j/sp_j)一定比(t_j/p_j)大。于是当前点左边的点就没有用了。因为直线斜率递增，所以最优解的位置递增。

于是只要用(Graham)算法维护凸壳即可。点排序后复杂度是(O(n))。

法3 分治dp

由于决策单调，计算(f[j],g[j])时还可以分治计算。

代码

//平衡树动态维护凸包
#include <bits/stdc++.h>
using namespace std;
#define mem(a,b) memset(a,b,sizeof(a))
#define rep(i,l,r) for(int i=l,_=r;i<_;++i)
#define per(i,l,r) for(int i=r-1,_=l;i>=_;--i)
typedef long long ll;
typedef double dd;
const int INF=0x3f3f3f3f;
const int N=50100;

struct nd{
	ll t,p;
	bool operator < (const nd&b)const{
		return t*b.p<p*b.t;
	}
}a[N];

struct DynmcCnvx{
	int rot,fa[N],c[N][2];
	dd lk[N],rk[N],x[N],y[N];
	void zigzag(int x,int &rot){
		int y=fa[x],z=fa[y];
		int p=(c[y][1]==x),q=p^1;
		if (y==rot) rot=x;
		else if (c[z][0]==y) c[z][0]=x; 
		else c[z][1]=x;
		fa[x]=z; fa[y]=x; fa[c[x][q]]=y;
		c[y][p]=c[x][q]; c[x][q]=y; 
	}
	void splay(int x,int &rot){
		while (x!=rot){
			int y=fa[x],z=fa[y];
			if (y!=rot) zigzag(((c[y][0]==x)^(c[z][0]==y))?x:y,rot);
			zigzag(x,rot);
		}
	}
	void insert(int &t,int anc,int now){//加入平衡树
		if (!t) t=now, fa[t]=anc;
		else insert(c[t][x[now]>x[t]],t,now);
	}
	void update(int t){//加入点(x[t],y[t])，维护上凸壳。
		splay(t,rot);
		if (c[t][0]){//向左求凸包
			int left=prev(rot);
			splay(left,c[rot][0]); c[left][1]=0;
			lk[t]=rk[left]=getk(left,t);
		}
		else lk[t]=INF;
		if (c[t][1]){//向右求凸包
			int right=succ(rot);
			splay(right,c[rot][1]); c[right][0]=0;
			rk[t]=lk[right]=getk(t,right);
		}
		else rk[t]=-INF;
		if (lk[t]<=rk[t]){//在原凸包内部的情况，直接删掉该点 
			rot=c[t][0]; c[rot][1]=c[t][1]; fa[c[t][1]]=rot; fa[rot]=0;
			lk[rot]=rk[c[t][1]]=getk(rot,c[t][1]);
		}
	}
	dd getk(int i,int j){//求斜率
		if (x[i]==x[j]) return -INF;
		return (y[j]-y[i])/(x[j]-x[i]);
	}
	int prev(int rot){//求可以和当前点组成凸包的右边第一个点
		int t=c[rot][0],tmp=t;
		while (t){
			if (getk(t,rot)<=lk[t]) tmp=t,t=c[t][1];
			else t=c[t][0];
		}
		return tmp;
	}
	int succ(int rot){//求可以和当前点组成凸包的左边第一个点
		int t=c[rot][1],tmp=t;
		while (t){
			if (getk(rot,t)>=rk[t]) tmp=t,t=c[t][0];
			else t=c[t][1];
		}
		return tmp;
	}
	int find(int t,dd k){//找到当前斜率的位置，即找到最优值
		if (!t) return 0;
		if (lk[t]>=k && k>=rk[t]) return t;
		return find(c[t][lk[t]>=k],k);
	}
  
	void Init(){
		rot=0;mem(fa,0);mem(c,0);
	}
	dd GetMax(dd a,dd b){//max{ax+by}
		int j=find(rot,-a/b);
		return a*x[j]+b*y[j];
	}
	void InsertPoint(int i,dd _x,dd _y){//插入点(x,y)
		x[i]=_x,y[i]=_y;
		insert(rot,0,i);
		update(i);
	}
}s;

dd st,sp;
dd f[N],g[N];
int ans[N],cnt;
int n;
int main(){
	ios::sync_with_stdio(false);
	cin.tie(0);
	cin>>n;
	rep(i,1,n+1)
		cin>>a[i].t>>a[i].p;
	sort(a+1,a+n+1);
	per(i,1,n+1){
		s.InsertPoint(i,-a[i].p,-a[i].t);
		f[i-1]=-s.GetMax(-(st+=a[i].t),sp+=a[i].p);
	}
	s.Init();
	rep(i,1,n){
		s.InsertPoint(i,a[i].p,a[i].t);
		g[i]=s.GetMax(-(st-=a[i].t),sp-=a[i].p);
		if(g[i]>f[i]) ans[cnt++]=i;
	}
	cout<<cnt<<endl;
	rep(i,0,cnt)cout<<ans[i]<<endl;
	return 0;
}

//直接维护凸包
#include <bits/stdc++.h>
using namespace std;
#define rep(i,l,r) for(int i=l,_=r;i<_;++i)
#define per(i,l,r) for(int i=r-1,_=l;i>=_;--i)
typedef long long ll;
const int N=50100;

struct nd{
	ll t,p;
	bool operator < (const nd&b)const{
		return t*b.p<p*b.t;
	}
}a[N];

struct Po{
	ll x,y;
	Po(ll x=0,ll y=0):x(x),y(y){}
	Po operator -(const Po&b)const {return Po(x-b.x,y-b.y);}
	Po operator +(const Po&b)const {return Po(x+b.x,y+b.y);}
	ll operator ^(const Po&b)const {return x*b.y-y*b.x;}
	ll operator *(const Po&b)const {return x*b.x+y*b.y;}
}p[N];
ll xmul(const Po&a,const Po&b,const Po&o){
	return (a-o)^(b-o);
}

struct DownCnvx{
	Po q[N];int top;
	//顺时针方向维护下凸壳
	void Insert(Po p){
		while(top && q[top].x<=p.x) --top;
		while(top>1 && xmul(q[top],p,q[top-1])>=0)--top;
		q[++top]=p;
	}
	ll GetMin(Po p){
		while(top>1 && p*q[top]>=p*q[top-1])--top;
		return q[top]*p;
	}
}d;

struct UpCnvx{
	Po q[N];int top;
	//顺时针方向维护上凸壳
	void Insert(Po p){
		while(top && q[top].x>=p.x) --top;
		while(top>1 && xmul(q[top],p,q[top-1])>=0)--top;
		q[++top]=p;
	}
	ll GetMax(Po p){
		while(top>1 && p*q[top]<=p*q[top-1])--top;
		return q[top]*p;
	}
}u;

ll st,sp;
ll f[N],g[N];
int cnt,ans[N];
int n;
int main(){
	ios::sync_with_stdio(false);
	cin.tie(0);
	cin>>n;
	rep(i,1,n+1)
		cin>>a[i].t>>a[i].p;
	sort(a+1,a+n+1);
	per(i,1,n+1){
		d.Insert(Po(a[i].p,a[i].t));
		f[i-1]=d.GetMin(Po(-(st+=a[i].t),sp+=a[i].p));
	}
	rep(i,1,n+1){
		u.Insert(Po(a[i].p,a[i].t));
		g[i]=u.GetMax(Po(-(st-=a[i].t),sp-=a[i].p));
		if(g[i]>f[i]) ans[cnt++]=i;
	}
	cout<<cnt<<endl;
	rep(i,0,cnt)cout<<ans[i]<<endl;
	return 0;
}

//分治优化
#include <bits/stdc++.h>
using namespace std;
#define rep(i,l,r) for(int i=l,_=r;i<_;++i)
#define per(i,l,r) for(int i=r-1,_=l;i>=_;--i)
typedef long long ll;
const ll LINF=0x3f3f3f3f3f3f3f3f;
const int N=50100;

struct nd{
	ll t,p;
	bool operator < (const nd&b)const{
		return t*b.p<p*b.t;
	}
}a[N];

ll st[N],sp[N];
ll f[N],g[N];
int cnt,ans[N];
int n;
void solveMax(int l,int r,int optL,int optR){
	if(l>r)return;
	int j=l+r>>1,u=optL;
	rep(i,optL,min(optR,j)+1){
		ll tmp=sp[j]*a[i].t-st[j]*a[i].p;
		if(tmp>g[j])g[j]=tmp,u=i;
	}
	solveMax(l, j-1, optL, u);
	solveMax(j+1, r, u, optR);
}
void solveMin(int l,int r,int optL,int optR){
	if(l>r)return;
	int j=l+r>>1,u=optL;
	rep(i,max(optL,j+1),optR+1){
		ll tmp=sp[j]*a[i].t-st[j]*a[i].p;
		if(tmp<f[j])f[j]=tmp,u=i;
	}
	solveMin(l, j-1, optL, u);
	solveMin(j+1, r, u, optR);
}
int main(){
	ios::sync_with_stdio(false);
	cin.tie(0);
	cin>>n;
	rep(i,1,n+1)
		cin>>a[i].t>>a[i].p;
	sort(a+1,a+n+1);
	per(i,1,n+1){
		st[i-1]=st[i]+a[i].t;sp[i-1]=sp[i]+a[i].p;
		f[i]=LINF;g[i]=-LINF;
	}
	solveMax(1,n,1,n);
	solveMin(1,n,1,n);
	rep(i,1,n)
		if(f[i]<g[i])ans[cnt++]=i;
	cout<<cnt<<endl;
	rep(i,0,cnt)cout<<ans[i]<<endl;
	return 0;
}