LOJ3059. 「HNOI2019」序列
给出序列\(a_i\),你要决定一个实数不下降序列\(b_i\),最小化\(\sum(a_i-b_i)^2\)。
有若干个询问,表示修改一个\(a_i\)之后答案。
\(n\le 10^5\)
参考题解:https://www.luogu.com.cn/blog/zhongyuwei/solution-p5294
考察对保序回归的进一步理解。
关于解的唯一性的证明:可以把\(\{a_i\},\{b_i\}\)分别看做高维向量,要最小化它们之间的距离。\(b_i\le b_{i+1}\)得到的空间是凸的,所以解唯一。
假设当前位置是\(x\)。考虑求出在搞完之后,单调栈中\(x\)所在区间的左右端点。求出区间\([1,x-1]\)和\([x+1,n]\)(从后往前)的单调栈。现在要把这两个单调栈和\(a_x\)合并。
一个算法是:枚举\(R_0\)表示后面单调栈的第\(R_0\)个,把\(x\)和\([R_0+1,tp_r]\)先合并,合并之后尝试去和前面的栈合并。如果合并出来了之后小于\(v(R_0)\),则合法。找到最大的\(R_0\),并且二分出最大的\(L_0\)(表示\([L_0+1,tp_l],x,[R_0+1,tp_r]\)合并起来大于\(v(L_0)\)),如果小于\(v(R_0)\),则得到了我们想要的东西,即\([L_0+1,tp_l],x,[R_0+1,tp_r]\)是最终\(x\)所在单调栈的区间。
接下来证明\(R_0\)的可二分性:也就是说如果\(R_0\)合法,\(R_0-1\)也合法。设两者求出的\(L_0\)分别为\(L,L'\)。
记\(v(l,r)\)表示\([l+1,tp_l],x,[r+1,tp_r]\)合并之后的值。\(L_0\)是最大的\(L_0\)满足\(v(L_0,R_0)>v(L_0)\),等价于使\(v(L_0,R_0)\)最大的\(L_0\)。
因为\(v(L',R_0)<v(L,R_0)<v(R_0)<v(R_0-1)\),\(v(L',R_0-1)\)是\(v(L',R_0),v(R_0)\)的带权平均数,所以\(v(L',R_0-1)<v(R_0-1)\)。
于是可以二分套二分,先二分\(R_0\)再二分\(L_0\)。时间\(O(n\lg^2 n)\)。
using namespace std;
#include <bits/stdc++.h>
#define N 100005
#define ll long long
#define ull unsigned long long
#define mo 998244353
ll sqr(ll x){return x*x%mo;}
ll qpow(ll x,ll y=mo-2){
ll r=1;
for (;y;y>>=1,x=x*x%mo)
if (y&1)
r=r*x%mo;
return r;
}
ll inv[N];
int n,m;
int a[N];
struct Query{int x,c,id;} q[N];
bool cmpq(Query a,Query b){return a.x<b.x;}
bool cmpf(ll s0,ll w0,ll s1,ll w1){return (ull)s0*(ull)w1<=(ull)s1*(ull)w0;}//用<=而不是<(可能会出现分子分母为0,此时直接返回true)
struct info{
ll w,v,v2,s,sw,sv,sv2;
} sl[N],sr[N];
int tpl,tpr;
info merge(info a,info b){
a.w+=b.w;
a.v+=b.v;
(a.v2+=b.v2)%=mo;
return a;
}
struct oper{
int t,ty;
info c;
} o[N*2];
int cnt;
int calcl(int R0,int c){
int l=0,r=tpl,_L0=0;
while (l<=r){
int L0=l+r>>1;
ll tmpv=c+sl[tpl].sv-sl[L0].sv+sr[tpr].sv-sr[R0].sv,tmpw=1+sl[tpl].sw-sl[L0].sw+sr[tpr].sw-sr[R0].sw;
if (cmpf(sl[L0].v,sl[L0].w,tmpv,tmpw)){
_L0=L0;
l=L0+1;
}
else
r=L0-1;
}
return _L0;
}
ll calcr(int i,int c){
int l=0,r=tpr,_R0=0,_L0=0;
while (l<=r){
int R0=l+r>>1,L0=calcl(R0,c);
ll tmpv=c+sl[tpl].sv-sl[L0].sv+sr[tpr].sv-sr[R0].sv,tmpw=1+sl[tpl].sw-sl[L0].sw+sr[tpr].sw-sr[R0].sw;
if (cmpf(tmpv,tmpw,sr[R0].v,sr[R0].w)){
_R0=R0,_L0=L0;
l=R0+1;
}
else
r=R0-1;
}
ll tmpv=c+sl[tpl].sv-sl[_L0].sv+sr[tpr].sv-sr[_R0].sv,tmpw=1+sl[tpl].sw-sl[_L0].sw+sr[tpr].sw-sr[_R0].sw;
ll tmpv2=sqr(c)+sl[tpl].sv2-sl[_L0].sv2+sr[tpr].sv2-sr[_R0].sv2;
ll res=(sl[_L0].s+sr[_R0].s-sqr(tmpv%mo)*inv[tmpw]%mo+tmpv2)%mo;
return res;
}
ll ans[N];
int main(){
freopen("in.txt","r",stdin);
scanf("%d%d",&n,&m);
for (int i=1;i<=n;++i)
scanf("%d",&a[i]);
for (int i=1;i<=m;++i)
scanf("%d%d",&q[i].x,&q[i].c),q[i].id=i;
inv[1]=1;
for (int i=2;i<=n;++i)
inv[i]=(mo-mo/i)*inv[mo%i]%mo;
for (int i=n;i>=1;--i){
info tmp={1,a[i],sqr(a[i])};
for (;tpr && cmpf(sr[tpr].v,sr[tpr].w,tmp.v,tmp.w);--tpr){
tmp=merge(tmp,sr[tpr]);
o[++cnt]={i,-1,sr[tpr]};
}
tmp.s=(sr[tpr].s-sqr(tmp.v%mo)*inv[tmp.w]+tmp.v2)%mo;
tmp.sw=sr[tpr].sw+tmp.w;
tmp.sv=sr[tpr].sv+tmp.v;
tmp.sv2=(sr[tpr].sv2+tmp.v2)%mo;
sr[++tpr]=tmp;
o[++cnt]={i,1};
}
ans[0]=sr[tpr].s;
sort(q+1,q+m+1,cmpq);
for (int i=1,j=1;i<=n;++i){
for (;cnt && o[cnt].t<=i;--cnt)
if (o[cnt].ty==-1)
sr[++tpr]=o[cnt].c;
else
--tpr;
for (;j<=m && q[j].x==i;++j)
ans[q[j].id]=calcr(i,q[j].c);
info tmp={1,a[i],sqr(a[i])};
for (;tpl && cmpf(tmp.v,tmp.w,sl[tpl].v,sl[tpl].w);--tpl)
tmp=merge(tmp,sl[tpl]);
tmp.s=(sl[tpl].s-sqr(tmp.v%mo)*inv[tmp.w]+tmp.v2)%mo;
tmp.sw=sl[tpl].sw+tmp.w;
tmp.sv=sl[tpl].sv+tmp.v;
tmp.sv2=(sl[tpl].sv2+tmp.v2)%mo;
sl[++tpl]=tmp;
}
for (int i=0;i<=m;++i)
printf("%lld\n",(ans[i]+mo)%mo);
return 0;
}