给定一个不降整数序列 \(1\le x_1\le x_2\le \cdots\le x_n\le q\)，请构造一个实数序列 \(y\) 满足 \(y_i\in [1,q]\)，\(y_i-y_{i-1}\in[a,b]\)，且最小化 \(\sum (y_i-x_i)^2\)，保证有解。
利用凸函数性质维护导数
我们设 \(dp_i(u)\) 表示对于所有的合法的 \(u\)，\(y_i=u\) 时 \(\sum_{j\le i}(y_j-x_j)^2\) 的最小值。
然后有转移式，设 \(dp_i(u)\) 的最小值在 \(v\) 处取到，那么 \(dp_{i+1}(u)\) 就应该是 \(dp_{i+1}(u)=\min_{d\in[a,b]} dp_{i}(u-d)+(u-x_j)^2\)。
我们发现 \(dp_i(u)\) 是有点棘手的，但是我们转而研究 \(dp_i'(u)\)。首先，我们发现 \(dp_i'(u)\) 是单调的，也就是 \(dp_i(u)\) 是凸的。考虑数学归纳法，在 \(dp_1(u)\)，因为只是一个开口向上的二次函数，所以其导数必然是单调增的。
然后如果 \(dp_{i-1}(u)\) 满足条件，就很好了，因为凸函数的转移是显然的。设 \(dp_{i-1}(u)\) 的最小值取在 \(k_i\)（注意因为有定义域限制，所以此处导数不一定为 \(0\)）
所能转移到的函数段单调下降，
所能转移到的函数段包含顶点，
所能转移到的函数段单调上升，
这是很好证明的，在凸函数上这些都是显然的。然后我们就发现，如果是导数的话：
这就很好了！因为我们发现，新的函数就是把原来的函数从中间断开，然后左边往右平移 \(a\)，右边往右平移 \(b\)，中间用 \(0\) 填满。这是单调不降的。然后整体加上单调上升的 \(2u-2x_i\)，我们发现 \(dp_i(u)\) 也是单调上升的。我们就证明了任意的 \(dp_i(u)\) 都是单调上升的。
然后我们发现，我们每次只是平移，然后加上一次函数，所以导函数也一直是一次函数，这就很好分段维护了。考虑分段维护当前所有的 \(dp_i'(x)\)，每次只需要找到零点断开，左右分别平移，中间安排 \(0\)，然后整体加上 \(2u-2x_i\)。每次最多增加两个段，所以最终的段个数不会超过 \(2n\)，复杂度 \(O(n^2)\)，是可以通过的。
讲一下卡了我很久的一个问题，就是我们的函数可能不是连续的。我一开始认为它一定是连续的，所以在发现函数不连续且没有 \(0\) 点的时候我就一直在找转移的问题，但实际上，我们的原函数本来就不是光滑的，插入一次函数 \(0u+0\) 之后导函数也不是连续的。但是我们的 \(0\) 点可以正常找，因为 \(0\) 点的定义并不是导数为 \(0\) 的点，而是原函数取到最小值的点。这时候如果相邻两段导函数一个右端点取值是负的，一个左端点取值是正的，就可以直接用它们交接的地方做 \(0\) 点，端点同理。
最后我们以 \(y_n=k_n\) 为基础倒推 \(y\)，因为 \(k_i\) 是使得当前取值最小的位置，而 \(dp_i(u)\) 还是个凸函数，所以很显然，尽可能的让 \(y_{i-1}\) 在不违反 \(y_i\) 条件的情况下最接近 \(k_{i-1}\) 即可。然后直接用得到的 \(y\) 计算答案。
const int N=6000;
const ld limit=0.000000001;
ll n,x[N+5],q,a,b;
ld y[N+5];
struct slope{
ld k,b,l,r;
slope(ld K,ld B,ld L,ld R){
k=K,b=B,l=L,r=R;
}
slope(){
k=0,b=0,l=0,r=0;
}
ld zerop(){
return -b/k;
}
};
slope dp[3*N+5],f[3*N+5];
ld k[N+5];
int m,cnt;
int main(){
scanf("%lld%lld%lld%lld",&n,&q,&a,&b);
rep(i,1,n)scanf("%lld",&x[i]);
dp[m=1]=slope(2,-2*x[1],1,q);
rep(i,2,n){
cnt=0;
k[i]=dp[1].l;
rep(j,1,m){
ld c=dp[j].zerop();
if(dp[j].l<c+limit && c<dp[j].r+limit)k[i]=c;
if(j!=1&&dp[j-1].k*dp[j].l+dp[j-1].b+limit<0){
if(dp[j].k*dp[j].l+dp[j].b>limit)k[i]=dp[j].l;
}
}
if(k[i]>q+limit)k[i]=q;
rep(j,1,m){
if(dp[j].r<k[i]+limit){//r<=k
f[++cnt]=slope(dp[j].k,dp[j].b-a*dp[j].k,dp[j].l+a,dp[j].r+a);
}else if(dp[j].l>k[i]+limit){//l>k
f[++cnt]=slope(dp[j].k,dp[j].b-b*dp[j].k,dp[j].l+b,dp[j].r+b);
}else{//l<=k<r
if(dp[j].l+limit<k[i]){//l<k
f[++cnt]=slope(dp[j].k,dp[j].b-a*dp[j].k,dp[j].l+a,k[i]+a);
}
f[++cnt]=slope(0,0,k[i]+a,k[i]+b);
f[++cnt]=slope(dp[j].k,dp[j].b-b*dp[j].k,k[i]+b,dp[j].r+b);
}
}
m=cnt;
rep(j,1,m){
dp[j]=f[j];
dp[j].k=dp[j].k+2;
dp[j].b=dp[j].b-2*x[i];
}
}
ld ans=dp[1].l;
rep(j,1,m){
ld c=dp[j].zerop();
if(dp[j].l<c+limit && c<dp[j].r+limit)ans=c;
if(j!=1&&dp[j-1].k*dp[j].l+dp[j-1].b+limit<0)
if(dp[j].k*dp[j].l+dp[j].b>limit)ans=dp[j].l;
}
if(ans>q+limit)ans=q;
per(i,1,n){
y[i]=ans;
if(ans-b+limit>k[i])ans=ans-b;
else if(ans-a<k[i]+limit)ans=ans-a;
else ans=k[i];
}
ld sum=0;
rep(i,1,n){
sum=sum+(x[i]-y[i])*(x[i]-y[i]);
printf("%.8Lf ",y[i]);
}
printf("\n%.8Lf\n",sum);
return 0;
}
//Nyn-siyn-hert
二次函数规划
金老师讲的做法，妙妙。
这个算法的核心是这样的。假设当 \(y_i=dp_i\) 时存在最优化 \(\sum_{j\le i}(x_j-y_j)^2\) 的一个方案。我们通过调整原先的答案依次求出 \(dp_i\)，然后同样选取结果倒推。
我们需要实现一个功能，一个函数 ld solve(int i,ld l,ld r,ld A,ld B)，表示返回一个 \(v\)，使得对于前 \(i\) 个位置，存在一个解 \(\{y\}\) 满足 \(y_i=v\) 并且 \(\sum_{j<i}(y_j-x_j)^2+Av^2+Bv\) 在 \(v\in[l,r]\) 的限制条件下 被该方案最优化。
我们发现，想要求出 \(dp_i\)，其实就是做一个 \(l=1,r=q\) 的规划，这个规划要最优化的是 \(\sum_{j\le i}(y_j-x_j)^2\)，忽略常数项（最优解的取值和常数项无关，即二次函数顶点横坐标和常数无关）也就是 \(\sum_{j<i}(y_j-x_j)^2+y_i^2-2x_iy_i\)。即 \(A=1,B=-2x_i\)。
然后考虑如何求取规划。
首先证明一个引理：
- 假设存在一个前 \(i\) 位的子问题的最优解 \(\{y\}\)，并且 \(y_i+b<cur\)（\(cur\) 是 \(Ay_{i+1}^2+By_{i}\) 的最优解），那么一定存在一个前 \(i+1\) 位的子问题的最优解 \(\{z\}\) 满足 \(z_i+b=z_{i+1}\)。
这个结论看起来是显然的，不过我们还是证一下。（金老师上课没讲，所以这个做法是我自己搞的，需要上一个做法的结论。不知道金老师的证法需不需要）
首先考虑最优解的 \(y_{i+1}\)，如果有多个解，我们选取 \(y_{i+1}\) 最大的一个，都相同就选取 \(y_{i}\) 最小的一个。
然后假设 \(y_{i}>y_{i+1}-b\)，那么因为第一个做法中转移函数的凸性，所以 \(y_i=y_{i+1}-b\) 一定不会更劣（其实是严格更优）。但是这就和当前解是 \(y_i\) 最小的最优解矛盾，所以得证。
同理，又有
- 假设存在一个前 \(i\) 位的子问题的最优解 \(\{y\}\)，\(y_i+a>cur\)（\(cur\) 是 \(Ay_{i+1}^2+By_{i}\) 的最优解），那么一定存在一个前 \(i+1\) 位的子问题的最优解 \(\{z\}\) 满足 \(z_i+a=z_{i+1}\)。
这个证明过程和上面是一样的，就不证了。
那么，现在我们尝试规划 \(y_i\)，先把前面的东西当作常数，单独找到 \(Ay_i^2+By_i\) 的最优解 \(y_i=cur\)。这是二次函数顶点，注意有边界。
然后，如果 \(cur\in[dp_{i-1}+a,dp_{i-1}+b]\)，那么我们就找到当前规划的最优解（前半段和后半段同时最优化），直接返回 \(cur\)。
如果 \(cur>dp_{i-1}+b\)，由引理可知一定存在一个解使得 \(y_{i-1}+b=y_{i}\)，那么就直接令 \(y_i=y_{i-1}+b\)，这样 \(\sum_{j<i}(y_j-x_j)^2+Ay_i^2+By_i\) 也就变成了 \(\sum_{j<i}(y_j-x_j)^2+A(y_{i-1}+b)^2+B(y_{i-1}+b)\)，忽略常数项即为 \(\sum_{j<i-1}(y_j-x_j)^2+(A+1)y_{i-1}^2+(B-2x_{i-1}+2bA)y_{i-1}\)。同时为了满足当前规划条件，还要 \(y_{i-1}\in[l-b,r-b]\)。这就是一个 \(i'=i-1\) 的子规划问题。递归求解直至和上一位不冲突或没有上一位即可。
每次最多倒退 \(O(n)\) 次，所以复杂度是 \(O(n^2)\)。
const ld eps=1e-8;
int n,x[6005],q,a,b;
ld dp[6005],y[6005],ans;
inline bool lt(ld x,ld y){
return x+eps<y;
}
inline bool gt(ld x,ld y){
return x>y+eps;
}
inline bool le(ld x,ld y){
return x<y+eps;
}
inline ld solve(int i,ld l,ld r,ld A,ld B){
if(lt(l,1))l=1;
if(lt(q,r))r=q;
ld cyr=-B/(2*A);
if(lt(r,cyr))cyr=r;
if(lt(cyr,l))cyr=l;
if(i==1||(le(dp[i-1]+a,cyr)&&le(cyr,dp[i-1]+b))){
return cyr;
}else if(gt(dp[i-1]+a,cyr)){
return solve(i-1,l-a,r-a,A+1,B-2*x[i-1]+2*a*A)+a;
}else{
return solve(i-1,l-b,r-b,A+1,B-2*x[i-1]+2*b*A)+b;
}
}
signed main(){
ios::sync_with_stdio(false);
cin.tie(0);
cin>>n>>q>>a>>b;
rp(i,n)cin>>x[i];
dp[1]=x[1];
rep(i,2,n){
if(le(dp[i-1]+a,x[i])&&le(x[i],dp[i-1]+b)){
dp[i]=x[i];
}else if(gt(dp[i-1]+a,x[i])){
dp[i]=solve(i,1,q,1,-2*x[i]);
}else{
dp[i]=solve(i,1,q,1,-2*x[i]);
}
}ld res=dp[n];
cout<<fixed<<setprecision(8);
per(i,1,n){
y[i]=res;
ans+=(res-x[i])*(res-x[i]);
if(gt(dp[i-1]+a,res))res=res-a;
else if(gt(res,dp[i-1]+b))res=res-b;
else res=dp[i-1];
}
rep(i,1,n)cout<<y[i]<<" ";
cout<<endl;
cout<<ans<<endl;
return 0;
}
//Nyn-siyn-hert
\(\text{FHQTreap}\) 优化导数维护
我们考虑使用类似 \(\text{FHQTreap}\) 的数据结构维护导数。数据结构是为算法服务的，我们可以调整数据结构从而更加契合我们的算法。
首先看到 split，这个在无旋平衡树很多时候是核心操作，但是这里我们只需要插入，还不是按照下标插入。而我们发现我们还有一个功能是找零点，那么为什么不把这两个功能合并起来呢？
具体而言就是，我们用两个函数 splitl 和 splitr，一个把右端点取值大于等于 \(0\) 的 \(split\) 到右边，这样先 splitr 一次得到 \(l\) 和 \(r\)，右边的子树右端点取值就都大于等于 \(0\)。然后 splitl 把左端点取值大于等于 \(0\) 的划分到 右边，这样我们再把 \(r\) \(split\) 成 \(tmp\) 和 \(r\)，我们会惊奇的发现：
如果 \(tmp\) 不为零，它一定是跨越正负的。也就是零点一定在 \(tmp\) 里面。
如果 \(tmp\) 为 \(0\)，那么意味着没有 \(0\) 点，那么 \(0\) 点就在左树和右树的交界点，这正好是右子树的最左点的左边界。我们就很好的把“找到目标位置”和“插入”直接合并起来了。
然后考虑我们需要进行的操作，一个是子树整体平移，一个是子树整体加一次函数。我们可以设计三个 \(tag\)，\(addk,addb,move\)。\(move\) 是把定义域向右平移 \(move\)，然后给 \(b\) 减去 \(kmove\)。\(add\) 就是简单的加法。这样就不可避免的涉及到下传多 \(tag\) 的顺序问题。这里，我们钦定先执行 \(move\)，再执行加法。
然后考虑合并 \(tag\)，原先有的 \(tag\) 是 \(add\) 和 \(mov\)，新增的是 \(add'\) 和 \(mov'\)。本来的顺序应该是 \(mov\rightarrow add\rightarrow mov'\rightarrow add'\)。
但是合并 \(tag\) 使得顺序变成了 \(mov\rightarrow mov'\rightarrow add\rightarrow add'\)。这样就相当于交换了 \(mov'\) 和 \(add\) 的位置。\(mov'\) 对 \(add\) 是没有贡献的，但是 \(add\) 对 \(mov'\) 是有贡献的。
具体而言，本来的 \(b\) 应该是 \(b+addb-mov'(k+addk)+addb'\)，但是现在却是 \(b+addb-mov'k+addb'\)，这就少减了一个 \(mov'addk\)，不过很好的是 \(add\) 对 \(mov\) 都是没有贡献的。所以直接把这个减在 \(addb\) 上就能消除交换顺序造成的影响了。
而我们只需要在 split 之后给 \(l\) 和 \(r\) 分别打上右移 \(a\) 和 \(b\) 的 \(tag\)。
然后考虑 merge 的问题。我们发现，我们的树还有一个好处，就是我们不需要上传什么东西。我们不需要维护区间信息。那么在没有 pushup 的情况下，我们的 merge 就会好些很多。
最后是插入新段。对于 \(tmp\) 为 \(0\) 的情况，可以直接把 \(tmp\) 作为新点插入。否则，我们可以开两个新点作为 \(tmp\) 的左右儿子，一个记录 \(tmp\) 零点左边的部分右移 \(a\) 位的结果，另一个记录右边部分右移 \(b\) 位的结果。然后就可以直接 merge，避免了开一堆新的单点子树然后套好多层 merge 的情况。
我们直接 split 出最终的零点，倒推答案。
const ld eps=1e-8;
int n,X[6005],q,a,b;
ld dp[6005],y[6005],ans;
inline bool lt(ld x,ld y){
return x+eps<y;
}
mt19937 rng(time(0));
struct node{
ld k,b,l,r,addk,addb;
int ls,rs,mov;
node(){}
node(ld _k,ld _b,ld _l,ld _r){
mov=0,addk=0,addb=0,ls=0,rs=0;
k=_k,b=_b,l=_l,r=_r;
}
}tr[114514];
int cnt=0,root;
#define ls(x) tr[x].ls
#define rs(x) tr[x].rs
inline void addtg(int x,ld k,ld b,int m){
if(!x)return;
tr[x].addb+=-m*tr[x].addk+b;
tr[x].addk+=k;
tr[x].mov+=m;
}
inline void pushdown(int x){
if(!x)return;
tr[x].l+=tr[x].mov;
tr[x].r+=tr[x].mov;
tr[x].b-=tr[x].mov*tr[x].k;
tr[x].k+=tr[x].addk;
tr[x].b+=tr[x].addb;
addtg(ls(x),tr[x].addk,tr[x].addb,tr[x].mov);
addtg(rs(x),tr[x].addk,tr[x].addb,tr[x].mov);
tr[x].addk=0,tr[x].addb=0,tr[x].mov=0;
}
inline void splitl(int x,int &l,int &r){
if(!x)return (void)(l=r=0);
pushdown(x);
ld val=tr[x].k*tr[x].l+tr[x].b;
if(le(0,val)){
r=x,splitl(ls(x),l,ls(x));
}else l=x,splitl(rs(x),rs(x),r);
}
inline void splitr(int x,int &l,int &r){
if(!x)return (void)(l=r=0);
pushdown(x);
ld val=tr[x].k*tr[x].r+tr[x].b;
if(lt(0,val)){
r=x,splitr(ls(x),l,ls(x));
}else l=x,splitr(rs(x),rs(x),r);
}
inline int merge(int l,int r){
pushdown(l);
pushdown(r);
if(!l)return r;
if(!r)return l;
int key=rng()%2;
if(key){
tr[r].ls=merge(l,tr[r].ls);
return r;
}else {
tr[l].rs=merge(tr[l].rs,r);
return l;
}
}
inline int find_first(int x){
return tr[x].ls?find_first(tr[x].ls):x;
}
signed main(){
ios::sync_with_stdio(false);
cin.tie(0);
cin>>n>>q>>a>>b;
rp(i,n)cin>>X[i];
cnt++,root=1,tr[0].l=q;
tr[1]=node(2,-2*X[1],1,q);
rep(i,2,n){
int l,x,r,R;
splitr(root,l,r),splitl(r,x,r);
addtg(l,0,0,a),addtg(r,0,0,b);
R=find_first(r);
if(!x){
x=++cnt,dp[i-1]=tr[R].l;
tr[x]=node(0,0,tr[R].l+a,tr[R].l+b);
root=merge(merge(l,x),r);
}else{
ld cyr=-tr[x].b/tr[x].k;
dp[i-1]=cyr;
if(lt(tr[x].l,cyr)){
tr[++cnt]=node(tr[x].k,tr[x].b-a*tr[x].k,tr[x].l+a,cyr+a);
tr[x].ls=cnt;
}
if(lt(cyr,tr[x].r)){
tr[++cnt]=node(tr[x].k,tr[x].b-b*tr[x].k,cyr+b,tr[x].r+b);
tr[x].rs=cnt;
}tr[x].k=0,tr[x].b=0,tr[x].l=cyr+a,tr[x].r=cyr+b;
root=merge(merge(l,x),r);
}
addtg(root,2,-2*X[i],0);
}
cout<<fixed<<setprecision(8);
int l,x,r,R;
splitr(root,l,r);
splitl(r,x,r);
R=find_first(r);
ld res=(x?(-tr[x].b/tr[x].k):tr[R].l);
if(lt(q,res))res=q;
per(i,1,n){
y[i]=res;
ans+=(res-X[i])*(res-X[i]);
if(gt(dp[i-1]+a,res))res=res-a;
else if(gt(res,dp[i-1]+b))res=res-b;
else res=dp[i-1];
}
rep(i,1,n)cout<<y[i]<<" ";
cout<<endl;
cout<<ans<<endl;
return 0;
}
//Nyn-siyn-hert
两个做法之间的共同点
实际上，我们发现，这道题最关键的地方就是“对前 \(i\) 个子问题分别求解”。我们每次先把前缀 \(i\) 当成独立的子问题求解，然后推出 \(i+1\)。
三个做法中，这个东西都至关重要。
在第一个和第三个做法中，这个点是原函数的顶点，导函数的零点，有了它才能决定新函数在哪里断开。
在第二个做法中，这个点是预处理的前缀问题的最优解。当一个前缀问题中，前缀的前缀前缀的最优解和前缀的后缀的最优规划互不冲突，就找到了当前子问题的解。然后就可以往后推了。
这是数学归纳思想，或者说递推思想在 \(\text{OI}\) 中重要的体现。
这道题我学到了什么
其实，最主要的是第三个做法，搞清楚了多个互相影响的 \(tag\) 同时下传的时候如何分析它们之间的影响并将其消除。在其他地方可谓用途很大。
然后，就是用分段函数 dp 的方法。不知道是不是我做题太少，没有遇到其他这样的题。
最后，是子问题规划的思想。但是在现在的我看来这个做法就像一个纸飞机，只有别人当面给我折叠一遍，我才能理解这个方法。否则这种巧妙的做法是我完全无法想到，下次遇到类似题目也难以再一次做出的。所以它对现在的我来说反而价值不大，不知道未来的自己能不能真正掌握这种思想。