APIO2010 特别行动队

此题是很好的凸壳优化的动态规划题。



【题意简述】：

给定一个含N个正数的序列，要求把序列分成若干个连续段。给每个连续段定义一个函数f，其值为A*sum^2 + B*sum + C。其中sum表示这个连续段的各个数字之和，A、B、C是三个给定的整数参数，其中A<0。要求这些连续段的f函数值之和最大。

【朴素动态规划】

用状态F[i]表示将序列的前i个数分成若干段可以得到的最大的f函数之和。则可得状态转移方程：

F[i] = max{F[j]+A*(S[i]-S[j])^2 + B*(S[i]-S[j]) + C}，0<=j<i。其中S[x]表示序列中第一个数到第x个数的和。

初始状态F[0]=0。

【凸壳优化的动态规划】

在求解状态F[i]时，对于任意两个决策x,y(0<=x<y<i),决策x优于决策y的充分必要条件是

F[x] + A*(S[i]-S[x])^2 + B*(S[i]-S[x]) + C > F[y] + A*(S[i]-S[y])^2 + B*(S[i]-S[y]) + C

变形得：

(F[y] + A*S[y]^2 – B*S[y]) – (F[x] + A*S[x]^2 – B*S[x]) < 2*A*S[i]*(S[y]-S[x])

令f(x)=F[x] + A*S[x]^2 – B*S[x]

则上式变为

f(y)-f(x) < 2*A*S[i]*(S[y]-S[x])

因为x<y且序列中都是正整数，所以有S[x]<S[y]，将S[y]-S[x]移项得

[f(y)-f(x)] / (S[y]-S[x]) < 2*A*S[i]

将{s[i],f(i)}看成平面上的点，则上式左边就等于{s[x],f(x)},{s[y],f(y)}所在直线的斜率。所以上式可以转化为：

K(x,y) < 2*A*S[i]

设F[i]的最优决策为k，对于F[i]的其他任意决策j：

当j<k,有K(j,k) >= 2*A*S[i];

当k<j,有K(k,j) < 2*A*S[i]。

所以，找最优决策的过程转化为在1..i-1号点中找到一点k，使得k和它左边点的连线斜率都大于2*A*S[i]，和他右边的点的连线斜率都小于等于2*A*S[i]。显然，点k必然在点1..i-1形成的上凸的凸壳上。

因此，我们用一个单调队列维护一个斜率递减的凸壳。每插入一个点，判断斜率是否递减，如果不递减则在队尾不断删除元素；否则将该点加到队尾。找最优决策时则判断队首两个点所在直线的斜率是否大于等于2*A*S[i]，如果是，则删除队首元素；否则返回队首元素所对应的决策。注意无论是插入还是查找最优决策时，如果队列中只剩一个点了，就直接插入或返回决策值就可以了。之所以可以这么做，是因为2*A*S[i]递减，所以最优决策值递增。

【程序】

#include<cstdio>#include<cstring>using namespace std;typedef long long big;const int MAXN=1000200,INF=0;int N,A,B,C; big S[MAXN],F[MAXN];inline double gx(int x){ return F[x]+A*S[x]*S[x]-B*S[x]; }inline big fx(big x){ return A*x*x + B*x + C; }struct queue{int head,tail;inline queue() { head=1; tail=0; V[head].k=INF; }struct node { int id; big x,y; double k; }V[MAXN];inline double getk(int p,int q){ return (gx(V[q].id)-gx(V[p].id))/(S[V[q].id]-S[V[p].id]); }void ins(int id,big x,big y){V[0].id=id; V[0].x=x; V[0].y=y;while (head<tail && V[tail].k <= getk(tail,0)) tail--;V[0].k=getk(tail,0); V[++tail]=V[0];}int get(double k){while (head<tail && V[head+1].k >= k) head++;return V[head].id;}}Q;void init(){scanf("%d/n%d%d%d/n",&N,&A,&B,&C);S[0]=0;for (int i=1; i<=N; i++){scanf("%lld",&S[i]);S[i]+=S[i-1];}scanf("/n");}void solve(){F[0]=0; Q.ins(0,S[0],F[0]);for (int i=1; i<=N; i++){int j=Q.get(2*A*S[i]); F[i]=F[j]+fx(S[i]-S[j]);Q.ins(i,S[i],F[i]);}}void print(){printf("%lld/n",F[N]);}int main(){freopen("commando.in","r",stdin); freopen("commando.out","w",stdout);init();solve();print();return 0;}

【评测结果】

开始测试 Ambusher

正在编译 Ambusher.commando ...成功

正在测试 Ambusher.commando.1(commando1.in)... 正确 0.003秒得分: 10

正在测试 Ambusher.commando.2(commando2.in)... 正确 0.004秒得分: 10

正在测试 Ambusher.commando.3(commando3.in)... 正确 0.004秒得分: 10

正在测试 Ambusher.commando.4(commando4.in)... 正确 0.006秒得分: 10

正在测试 Ambusher.commando.5(commando5.in)... 正确 0.006秒得分: 10

正在测试 Ambusher.commando.6(commando6.in)... 正确 0.047秒得分: 10

正在测试 Ambusher.commando.7(commando7.in)... 正确 0.222秒得分: 10

正在测试 Ambusher.commando.8(commando8.in)... 正确 0.372秒得分: 10

正在测试 Ambusher.commando.9(commando9.in)... 正确 0.432秒得分: 10

正在测试 Ambusher.commando.10(commando10.in)... 正确 0.453秒得分: 10

Ambusher.commando 的总分: 100

Ambusher 的总分: 100