CF280E - Sequence Transformation

给定一个不降整数序列 $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$，保证有解。

利用凸函数性质维护导数

我们设 $d{p}_i(u)$ 表示对于所有的合法的 $u$，$y_i=u$ 时 $\sum _{j\le i}(y_j-x_j{)}^2$ 的最小值。

然后有转移式，设 $d{p}_i(u)$ 的最小值在 $v$ 处取到，那么 $d{p}_{i+1}(u)$ 就应该是 $d{p}_{i+1}(u)=\underset{d\in [a,b]}{min}d{p}_i(u-d)+(u-x_j{)}^2$。

我们发现 $d{p}_i(u)$ 是有点棘手的，但是我们转而研究 $d{p}_i^{\prime }(u)$。首先，我们发现 $d{p}_i^{\prime }(u)$ 是单调的，也就是 $d{p}_i(u)$ 是凸的。考虑数学归纳法，在 $d{p}_1(u)$，因为只是一个开口向上的二次函数，所以其导数必然是单调增的。

然后如果 $d{p}_{i-1}(u)$ 满足条件，就很好了，因为凸函数的转移是显然的。设 $d{p}_{i-1}(u)$ 的最小值取在 $k_i$（注意因为有定义域限制，所以此处导数不一定为 $0$）

所能转移到的函数段单调下降，

$$d{p}_i(u)=d{p}_{i-1}(u-a)+(u-x_i{)}^2(u<k_i+a)$$

所能转移到的函数段包含顶点，

$$d{p}_i(u)=d{p}_{i-1}(k_i)+(u-x_i{)}^2(k_i+a\le u\le k_i+b)$$

所能转移到的函数段单调上升，

$$d{p}_i(u)=d{p}_{i-1}(u-b)+(u-x_i{)}^2(u>k_i+b)$$

这是很好证明的，在凸函数上这些都是显然的。然后我们就发现，如果是导数的话：

$$d{p}_i^{\prime }(u)=d{p}_{i-1}^{\prime }(u-a)+2u-2x_i(u<k_i+a)$$

$$d{p}_i^{\prime }(u)=2u-2x_i(k_i+a\le u\le k_i+b)$$

$$d{p}_i^{\prime }(u)=d{p}_{i-1}^{\prime }(u-b)+2u-2x_i(u>k_i+b)$$

这就很好了！因为我们发现，新的函数就是把原来的函数从中间断开，然后左边往右平移 $a$，右边往右平移 $b$，中间用 $0$ 填满。这是单调不降的。然后整体加上单调上升的 $2u-2x_i$，我们发现 $d{p}_i(u)$ 也是单调上升的。我们就证明了任意的 $d{p}_i(u)$ 都是单调上升的。

然后我们发现，我们每次只是平移，然后加上一次函数，所以导函数也一直是一次函数，这就很好分段维护了。考虑分段维护当前所有的 $d{p}_i^{\prime }(x)$，每次只需要找到零点断开，左右分别平移，中间安排 $0$，然后整体加上 $2u-2x_i$。每次最多增加两个段，所以最终的段个数不会超过 $2n$，复杂度 $O(n^2)$，是可以通过的。

讲一下卡了我很久的一个问题，就是我们的函数可能不是连续的。我一开始认为它一定是连续的，所以在发现函数不连续且没有 $0$ 点的时候我就一直在找转移的问题，但实际上，我们的原函数本来就不是光滑的，插入一次函数 $0u+0$ 之后导函数也不是连续的。但是我们的 $0$ 点可以正常找，因为 $0$ 点的定义并不是导数为 $0$ 的点，而是原函数取到最小值的点。这时候如果相邻两段导函数一个右端点取值是负的，一个左端点取值是正的，就可以直接用它们交接的地方做 $0$ 点，端点同理。

最后我们以 $y_n=k_n$ 为基础倒推 $y$，因为 $k_i$ 是使得当前取值最小的位置，而 $d{p}_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=d{p}_i$ 时存在最优化 $\sum _{j\le i}(x_j-y_j{)}^2$ 的一个方案。我们通过调整原先的答案依次求出 $d{p}_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+A{v}^2+Bv$ 在 $v\in [l,r]$ 的限制条件下 被该方案最优化。

我们发现，想要求出 $d{p}_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_i{y}_i$。即 $A=1,B=-2x_i$。

然后考虑如何求取规划。

首先证明一个引理：

• 假设存在一个前 $i$ 位的子问题的最优解 $\{y\}$，并且 $y_i+b<cur$（$cur$ 是 $A{y}_{i+1}^2+B{y}_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$ 是 $A{y}_{i+1}^2+B{y}_i$ 的最优解），那么一定存在一个前 $i+1$ 位的子问题的最优解 $\{z\}$ 满足 $z_i+a=z_{i+1}$。

这个证明过程和上面是一样的，就不证了。

那么，现在我们尝试规划 $y_i$，先把前面的东西当作常数，单独找到 $A{y}_i^2+B{y}_i$ 的最优解 $y_i=cur$。这是二次函数顶点，注意有边界。

然后，如果 $cur\in [d{p}_{i-1}+a,d{p}_{i-1}+b]$，那么我们就找到当前规划的最优解（前半段和后半段同时最优化），直接返回 $cur$。

如果 $cur>d{p}_{i-1}+b$，由引理可知一定存在一个解使得 $y_{i-1}+b=y_i$，那么就直接令 $y_i=y_{i-1}+b$，这样 $\sum _{j<i}(y_j-x_j{)}^2+A{y}_i^2+B{y}_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^{\prime }=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$，新增的是 $ad{d}^{\prime }$ 和 $mo{v}^{\prime }$。本来的顺序应该是 $mov\to add\to mo{v}^{\prime }\to ad{d}^{\prime }$。

但是合并 $tag$ 使得顺序变成了 $mov\to mo{v}^{\prime }\to add\to ad{d}^{\prime }$。这样就相当于交换了 $mo{v}^{\prime }$ 和 $add$ 的位置。$mo{v}^{\prime }$ 对 $add$ 是没有贡献的，但是 $add$ 对 $mo{v}^{\prime }$ 是有贡献的。

具体而言，本来的 $b$ 应该是 $b+addb-mo{v}^{\prime }(k+addk)+add{b}^{\prime }$，但是现在却是 $b+addb-mo{v}^{\prime }k+add{b}^{\prime }$，这就少减了一个 $mo{v}^{\prime }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 的方法。不知道是不是我做题太少，没有遇到其他这样的题。

最后，是子问题规划的思想。但是在现在的我看来这个做法就像一个纸飞机，只有别人当面给我折叠一遍，我才能理解这个方法。否则这种巧妙的做法是我完全无法想到，下次遇到类似题目也难以再一次做出的。所以它对现在的我来说反而价值不大，不知道未来的自己能不能真正掌握这种思想。