BZOJ1492：[NOI2007]货币兑换Cash （CDQ分治+斜率优化DP/平衡树维护凸壳）

题目传送门：http://www.lydsy.com/JudgeOnline/problem.php?id=1492

题目分析：被一道水题坑了两天……

首先朴素的O(n2)的DP是很好想的：记f[i]表示第i天不持有任何金券所能获得的最多人民币，很明显f[i]可以由f[i-1]更新而来（就是在1~i-1天就已经卖完所有金券，然后第i天不进行任何操作），还可以枚举一个j，看一下在第j天买入金券，在第i天的时候卖出金券能否获得更多收益。由于每一次操作都必定花光所有的人民币或卖光所有金券，于是：

f[i]=max(f[i],f[j]Rate[j]Rate[j]A[j]+B[j]A[i]+f[j]Rate[j]A[j]+B[j]B[i])

接下来我们要找到原版的题目，然后看一下n的范围：n<=105（BZOJ上没有n的范围），于是要考虑优化。设：

P[j]=f[j]Rate[j]Rate[j]A[j]+B[j],Q[j]=−f[j]Rate[j]A[j]+B[j]

（注意Q[j]右边的那个地方有个负号）

这样我们只要算出了f[j]，就可以算出P[j],Q[j]，而我们要求的就是P[j]A[i]−Q[j]B[i]的最大值，即(A[i]B[i]P[j]−Q[j])B[i]的最大值，然后这就是一个很明显的斜率优化。我们手推一下式子就可以看出要维护一个斜率单调上升的下凸壳，并在查找答案的时候找出一条最左边的有向线段j->k使得其斜率大于A[i]B[i]，这样j便是最优答案。

这说起来很简单，但我们发现P[j]并不一定单调递增，而且A[i]B[i]也不一定单调递增，于是就要用CDQ分治来处理。分治的时候先递归左边算出左半部分的f值，然后将左半部分按P值升序排序，构出凸壳；再将右半部分按A[i]B[i]的大小排序，从左往右扫一遍凸壳即可。至于对P的排序可以用归并排序，对A[i]B[i]的排序可以在主函数里先排一次序，再一边分治一边进行归并排序的逆过程。由于每一层的时间是线性的，所以总时间是O(nlog(n))的。

CODE：

#include<iostream>
#include<string>
#include<cstring>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<stdio.h>
#include<algorithm>
using namespace std;

const int maxn=100100;
const double eps=1e-7;

double f[maxn];
int que[maxn];
int head,tail;

int pid[maxn];
int cid[maxn];
int temp[maxn];
int cnt;

double A[maxn];
double B[maxn];
double R[maxn];

double P[maxn];
double Q[maxn];
int s,n;

bool Comp(int x,int y)
{
return A[x]*B[y]<A[y]*B[x];
}

void Push(int x)
{
que[++tail]=x;
while (head+1<tail)
{
int y=que[tail-1];
int z=que[tail-2];
double Left=(Q[x]-Q[y])*(P[y]-P[z]);
double Right=(Q[y]-Q[z])*(P[x]-P[y]);
if (Left-Right>eps) break;
tail--;
que[tail]=x;
}
}

bool Check(int k)
{
while (head<tail)
{
int x=que[head];
int y=que[head+1];
double Left=A[k]*(P[y]-P[x]);
double Right=B[k]*(Q[y]-Q[x]);
if (Left-Right<-eps) break;
head++;
}
int j=que[head];
f[k]=max(f[k],P[j]*A[k]-Q[j]*B[k]);
}

void CDQ(int Le,int Ri)
{
if (Le==Ri)
{
f[Le]=max(f[Le],f[Le-1]);
P[Le]=f[Le]*R[Le]/(R[Le]*A[Le]+B[Le]);
Q[Le]=-f[Le]/(R[Le]*A[Le]+B[Le]);
return;
}
int mid=(Le+Ri)>>1;

cnt=0;
for (int i=Le; i<=Ri; i++) if (cid[i]<=mid) temp[++cnt]=cid[i];
for (int i=Le; i<=Ri; i++) if (cid[i]>mid) temp[++cnt]=cid[i];
for (int i=1; i<=cnt; i++) cid[Le+i-1]=temp[i];
CDQ(Le,mid);

head=1,tail=0;
for (int i=Le; i<=mid; i++) Push(pid[i]);
for (int i=mid+1; i<=Ri; i++) Check(cid[i]);

CDQ(mid+1,Ri);
int Left=Le,Right=mid+1;
cnt=0;
while ( Left<=mid || Right<=Ri )
{
if (Left>mid)
{
temp[++cnt]=pid[Right++];
continue;
}
if (Right>Ri)
{
temp[++cnt]=pid[Left++];
continue;
}
int x=pid[Left],y=pid[Right];
if (P[x]<P[y]) temp[++cnt]=pid[Left++];
else temp[++cnt]=pid[Right++];
}
for (int i=1; i<=cnt; i++) pid[Le+i-1]=temp[i];
}

int main()
{
freopen("cash.in","r",stdin);
freopen("cash.out","w",stdout);

scanf("%d%d",&n,&s);
for (int i=1; i<=n; i++)
scanf("%lf%lf%lf",&A[i],&B[i],&R[i]),cid[i]=pid[i]=i;
sort(cid+1,cid+n+1,Comp);

f[0]=(double)s;
CDQ(1,n);
printf("%.3lf\n",f
);

return 0;
}

其实也可以根本不用CDQ分治，直接用平衡树维护不断加点之后的凸壳即可，举个例子：



假设现在凸壳上的点有A,B,C,D，我们用一棵treap按X为第一关键字，Y为第二关键字维护。尝试加进点F时，先找到它的前驱(B)和后继(C)，看一下折线BFC是否破坏凸性，如果不破坏的话就将F加进去，然后尝试不断弹出它前后的点。

虽然代码看上去很简单，但其实并不是这么好写（而且我一开始还写错了很多地方）。由于treap存的是点的坐标，而二叉查找的时候是要通过该点和上一个点的斜率来判断向左走还是向右走，于是我们在一边二叉走的时候还要一边维护前驱。如果某个点P有左儿子，那么它的前驱就是它左儿子中最右边的点，否则我们要找到一个最近的P的祖先F（或者P自身），使得F是它父亲的右儿子，此时F的父亲便是P的前驱。这个F可以在一边二叉查找的时候一边维护，如果F不存在，则P没有前驱。

我最初的代码犯了很多错，比如Q[j]右边的式子没有负号；想要维护某棵子树最右边的点的X和Y值，结果莫名写成了维护一棵子树中X,Y的最大值QAQ；二叉查找的时候，一个点没有前驱就不往它的右儿子走……之类的。改了之后交上去还是80，后来看讨论发现可能有重复的点，而我的treap虽然能处理X坐标相同，共线的情况，却处理不了这种情况。特判了一下就A了，改了我好久……

CODE：

#include<iostream>
#include<string>
#include<cstring>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<stdio.h>
#include<algorithm>
using namespace std;

const int maxn=100100;
const double eps=1e-8;

const long long M1=998244553;
const long long M2=1000000007;
const long long M3=1333333331;
typedef long long LL;

struct Tnode
{
double valX,valY,maxX,maxY;
int fix;
Tnode *lson,*rson;
void Up()
{
maxX=valX,maxY=valY;
if (rson) maxX=rson->maxX,maxY=rson->maxY;
}
} tree[maxn];
Tnode *Root=NULL;
int cur=-1;

double f[maxn];
double a[maxn];
double b[maxn];
double r[maxn];

LL seed;
int s,n;

int Rand()
{
seed=(seed*M1+M2)%M3;
return (int)seed;
}

Tnode *New_node(double x,double y)
{
cur++;
tree[cur].valX=tree[cur].maxX=x;
tree[cur].valY=tree[cur].maxY=y;
tree[cur].fix=Rand();
tree[cur].lson=tree[cur].rson=NULL;
return tree+cur;
}

double Abs(double x)
{
if (x>=0.0) return x;
return -x;
}

void Right_turn(Tnode *&P)
{
Tnode *W=P->lson;
P->lson=W->rson;
W->rson=P;
P=W;
P->rson->Up();
P->Up();
}

void Left_turn(Tnode *&P)
{
Tnode *W=P->rson;
P->rson=W->lson;
W->lson=P;
P=W;
P->lson->Up();
P->Up();
}

void Insert(Tnode *&P,double x,double y)
{
if (!P) P=New_node(x,y);
else
if ( x-P->valX<-eps || ( Abs(x-P->valX)<=eps && y-P->valY<-eps ) )
{
Insert(P->lson,x,y);
if ( P->lson->fix < P->fix ) Right_turn(P);
else P->Up();
}
else
{
Insert(P->rson,x,y);
if ( P->rson->fix < P->fix ) Left_turn(P);
else P->Up();
}
}

Tnode *Find(Tnode *P,Tnode *F,double c,Tnode *op)
{
if (!P) return op;
double lx,ly;
if (P->lson) lx=P->lson->maxX,ly=P->lson->maxY;
else
if (F) lx=F->valX,ly=F->valY;
else return Find(P->rson,P,c,P);
if ( (P->valY-ly)-c*(P->valX-lx)<-eps ) return Find(P->rson,P,c,P);
else return Find(P->lson,F,c,op);
}

bool Check(Tnode *P,double x,double y)
{
if (!P) return false;
if ( Abs(P->valX-x)<=eps && Abs(P->valY-y)<=eps ) return true;
if ( x-P->valX<-eps || ( Abs(x-P->valX)<=eps && y-P->valY<-eps ) )
return Check(P->lson,x,y);
return Check(P->rson,x,y);
}

Tnode *Get_prev(Tnode *P,double x,double y,Tnode *op)
{
if (!P) return op;
if ( P->valX-x<-eps || ( Abs(x-P->valX)<=eps && P->valY-y<-eps ) )
return Get_prev(P->rson,x,y,P);
return Get_prev(P->lson,x,y,op);
}

Tnode *Get_succ(Tnode *P,double x,double y,Tnode *op)
{
if (!P) return op;
if ( P->valX-x>eps || ( Abs(x-P->valX)<=eps && P->valY-y>eps ) )
return Get_succ(P->lson,x,y,P);
return Get_succ(P->rson,x,y,op);
}

bool Judge(double ax,double ay,double bx,double by,double cx,double cy)
{
return ( (cy-by)*(bx-ax)-(by-ay)*(cx-bx)>eps );
}

void Delete(Tnode *&P,double x,double y)
{
if ( Abs(P->valX-x)<=eps && Abs(P->valY-y)<=eps )
if (P->lson)
if (P->rson)
{
if ( P->lson->fix < P->rson->fix )
{
Right_turn(P);
Delete(P->rson,x,y);
}
else
{
Left_turn(P);
Delete(P->lson,x,y);
}
P->Up();
}
else P=P->lson;
else P=P->rson;
else
{
if ( x-P->valX<-eps || ( Abs(x-P->valX)<=eps && y-P->valY<-eps ) )
Delete(P->lson,x,y);
else Delete(P->rson,x,y);
P->Up();
}
}

int main()
{
freopen("cash.in","r",stdin);
freopen("cash.out","w",stdout);

scanf("%d%d",&n,&s);
seed=s;
for (int i=1; i<=n; i++) scanf("%lf%lf%lf",&a[i],&b[i],&r[i]);

f[1]=s;
double p=f[1]*r[1]/(r[1]*a[1]+b[1]);
double q=-f[1]/(r[1]*a[1]+b[1]);
Insert(Root,p,q);

for (int i=2; i<=n; i++)
{
Tnode *P=Find(Root,NULL,a[i]/b[i],NULL);
f[i]=max(f[i-1],P->valX*a[i]-P->valY*b[i]);
p=f[i]*r[i]/(r[i]*a[i]+b[i]);
q=-f[i]/(r[i]*a[i]+b[i]);
if ( Check(Root,p,q) ) continue;

bool flag=true;
Tnode *Prev=Get_prev(Root,p,q,NULL);
Tnode *Succ=Get_succ(Root,p,q,NULL);
if ( Prev && Succ ) flag=Judge(Prev->valX,Prev->valY,p,q,Succ->valX,Succ->valY);

if (flag)
{
Insert(Root,p,q);
Tnode *tp1,*tp2;
while (1)
{
tp1=Get_prev(Root,p,q,NULL);
if (!tp1) break;
tp2=Get_prev(Root,tp1->valX,tp1->valY,NULL);
if (!tp2) break;
if ( Judge(tp2->valX,tp2->valY,tp1->valX,tp1->valY,p,q) ) break;
Delete(Root,tp1->valX,tp1->valY);
}

while (1)
{
tp1=Get_succ(Root,p,q,NULL);
if (!tp1) break;
tp2=Get_succ(Root,tp1->valX,tp1->valY,NULL);
if (!tp2) break;
if ( Judge(p,q,tp1->valX,tp1->valY,tp2->valX,tp2->valY) ) break;
Delete(Root,tp1->valX,tp1->valY);
}
}
}

printf("%.3lf\n",f
);
return 0;
}