原题题解和数据下载 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_{i\in S}t_i} {\sum_{i\in S}p_i}>r_j
\]
\]
等价于\(\sum_{i\in S}t_i- {\sum_{i\in S}p_i} \cdot r_j>0\)。也就是：
\[\sum_{i\in S}(t_i-p_i \cdot r_j)>0 \Rightarrow
\sum_{i\in S}t_i\cdot sp_j-p_i\cdot st_j>0
\]
\sum_{i\in S}t_i\cdot sp_j-p_i\cdot st_j>0
\]
贪心地找最大的 n-j 个加起来，复杂度是\(O(n^2\log n)\)。
法1 平衡树维护动态凸包
一个优化的方法是，考虑\(z=t_i\cdot sp_j-p_i\cdot st_j\)在\(1..j\)里最大值\(g[j]\) 和\(j+1..n\)里的最小值\(f[j]\)。
如果\(f[j]\) 小于\(g[j]\)，那么意味着从\(j+1..n\)中取出最小值换为\(1..j\)中的最大值，可以更优。
\(g[j]=\max\{a_j\cdot x_i+b_j\cdot y_i\}\)实际上就是斜率优化的形式。
这里\(x_i=p_i,y_i=t_i,a_j=-st_j,b_j=sp_j\)。
直线为\(y=st_j/sp_j\cdot x+z/sp_j\)。
所以每个点坐标为\((p_i,t_i)\)，要找一个点，斜率为\(-a_j/b_j\)的直线经过它时，纵截距最大。需要维护一个上凸壳。x不是单调的，需要用平衡树来维护这个凸壳。
\(f[j]\)同理求。总复杂度\(O(n\log 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;
}