使用矩阵乘法+线段树解决区间历史和问题的一种通用解法

前言

一般解决普通的区间历史和，只需要定义辅助$c=hs-t\cdot a$，$hs$是历史和，$a$是区间和，$t$是时间戳，维护$a,c$数组的区间加即可

但是如果题目更复杂一点，就要设置不同的tag，考虑下放顺序和影响，非常费脑子

一种无脑的方法就是构造矩阵，如果单纯使用普通的矩阵乘法，那么总复杂度会多$C^3$，其中$C$为向量长度

本文的目的就是优化矩阵乘法的过程，实际上就是对矩阵乘法循环展开，只不过个人认为能更优雅一点罢了

P8868 [NOIP2022] 比赛

题意给出$a,b$两个数组，多次询问，$[L,R]$的所有子区间的$\max a_i\cdot \max b_j$，也就是求$$\sum _{L\le l\le r\le R}\underset{l\le i\le r}{\max }{a}_i\cdot \underset{l\le i\le r}{\max }{b}_i$$

考虑对所有询问离线，从左到右扫一遍，维护以$i$为右端点的答案，询问就是查询区间$\max a\cdot \max b$的历史和

对$a,b$分别维护两个单调栈即可用线段树更新

线段树每个节点维护$[a,b,ab,c,len]$表示区间$a_i$的和，$b_i$的和，$a_i{b}_i$的和，$a_i{b}_i$的历史和，区间长度

那么对于区间$a$加$k$，有
$${\left(\begin{array}{c}a\\ b\\ ab\\ c\\ len\end{array}\right)}^T\cdot \left(\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& k& 0& 0\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 1& 0\\ k& 0& 0& 0& 1\end{array}\right)={\left(\begin{array}{c}a+k\cdot len\\ b\\ ab+kb\\ c\\ len\end{array}\right)}^T$$
同理，区间$b$加$k$,有
$$\left(\begin{array}{ccccc}1& 0& k& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1\\ 0& k& 0& 0& 1\end{array}\right)$$
更新区间历史和
$$\left(\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1& 1& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right)$$
虽然矩阵乘法是$5^3$，但是可以发现很多状态是一直为$0$，也就是没有用的，要想找到这些状态，我们只需要将所有矩阵初始不为0的状态设为1，跑一遍传递闭包，就可以知道所有状态，这里25个状态一共只有14个状态合法(实际上可以9个，主对角线恒为1)

然后我们可以打表来代替手写，打表代码如下

void solve(){vector<vector<int>> f={{1,0,1,0,0},{0,1,1,0,0},{0,0,1,1,0},{0,0,0,1,0},{1,1,0,0,1}};const int n=f.size();for(int k=0;k<n;k++){for(int i=0;i<n;i++){for(int j=0;j<n;j++)f[i][j]|=f[i][k]&f[k][j];}}string s="int ";for(int i=0;i<n;i++){for(int j=0;j<n;j++){if(f[i][j])s+="x"+to_string(i)+to_string(j)+",";}}s.pop_back();s+=";\n";s+="\n\nMatrix * Matrix:\n";for(int i=0;i<n;i++){for(int j=0;j<n;j++){if(!f[i][j])continue;s+="res.x"+to_string(i)+to_string(j)+"=";for(int k=0;k<n;k++){if(f[i][k]&f[k][j])s+="a.x"+to_string(i)+to_string(k)+"*b.x"+to_string(k)+to_string(j)+"+";}s.pop_back();s+=";\n";}}s+="\n\nVec * Matrix:\n";for(int i=0;i<n;i++){s+="res.x0"+to_string(i)+"=";for(int j=0;j<n;j++){if(f[j][i]){s+="a.x0"+to_string(j)+"*b.x"+to_string(j)+to_string(i)+"+";}}s.pop_back();s+=";\n";}s.pop_back();cout<<s;
}
/*
打印结果：
int x00,x02,x03,x11,x12,x13,x22,x23,x33,x40,x41,x42,x43,x44;Matrix * Matrix:
res.x00=a.x00*b.x00;
res.x02=a.x00*b.x02+a.x02*b.x22;
res.x03=a.x00*b.x03+a.x02*b.x23+a.x03*b.x33;
res.x11=a.x11*b.x11;
res.x12=a.x11*b.x12+a.x12*b.x22;
res.x13=a.x11*b.x13+a.x12*b.x23+a.x13*b.x33;
res.x22=a.x22*b.x22;
res.x23=a.x22*b.x23+a.x23*b.x33;
res.x33=a.x33*b.x33;
res.x40=a.x40*b.x00+a.x44*b.x40;
res.x41=a.x41*b.x11+a.x44*b.x41;
res.x42=a.x40*b.x02+a.x41*b.x12+a.x42*b.x22+a.x44*b.x42;
res.x43=a.x40*b.x03+a.x41*b.x13+a.x42*b.x23+a.x43*b.x33+a.x44*b.x43;
res.x44=a.x44*b.x44;Vec * Matrix:
res.x00=a.x00*b.x00+a.x04*b.x40;
res.x01=a.x01*b.x11+a.x04*b.x41;
res.x02=a.x00*b.x02+a.x01*b.x12+a.x02*b.x22+a.x04*b.x42;
res.x03=a.x00*b.x03+a.x01*b.x13+a.x02*b.x23+a.x03*b.x33+a.x04*b.x43;
res.x04=a.x04*b.x44;
*/

然后我们就可以利用打印结果，快速套上线段树板子

struct Matrix{int x00,x02,x03,x11,x12,x13,x22,x23,x33,x40,x41,x42,x43,x44;
};
struct Vec{int x00,x01,x02,x03,x04;
};
Matrix operator*(const Matrix &a,const Matrix &b){Matrix res;res.x00=a.x00*b.x00;res.x02=a.x00*b.x02+a.x02*b.x22;res.x03=a.x00*b.x03+a.x02*b.x23+a.x03*b.x33;res.x11=a.x11*b.x11;res.x12=a.x11*b.x12+a.x12*b.x22;res.x13=a.x11*b.x13+a.x12*b.x23+a.x13*b.x33;res.x22=a.x22*b.x22;res.x23=a.x22*b.x23+a.x23*b.x33;res.x33=a.x33*b.x33;res.x40=a.x40*b.x00+a.x44*b.x40;res.x41=a.x41*b.x11+a.x44*b.x41;res.x42=a.x40*b.x02+a.x41*b.x12+a.x42*b.x22+a.x44*b.x42;res.x43=a.x40*b.x03+a.x41*b.x13+a.x42*b.x23+a.x43*b.x33+a.x44*b.x43;res.x44=a.x44*b.x44;return res;
}
Vec operator * (const Vec &a,const Matrix &b){Vec res;res.x00=a.x00*b.x00+a.x04*b.x40;res.x01=a.x01*b.x11+a.x04*b.x41;res.x02=a.x00*b.x02+a.x01*b.x12+a.x02*b.x22+a.x04*b.x42;res.x03=a.x00*b.x03+a.x01*b.x13+a.x02*b.x23+a.x03*b.x33+a.x04*b.x43;res.x04=a.x04*b.x44;return res;
}
Vec operator +(const Vec &a,const Vec &b){return {a.x00+b.x00,a.x01+b.x01,a.x02+b.x02,a.x03+b.x03,a.x04+b.x04};
}
Matrix I={1,0,0,1,0,0,1,0,1,0,0,0,0,1};
Matrix C={1,0,0,1,0,0,1,1,1,0,0,0,0,1};
Matrix getA(int k){return Matrix{1,0,0,1,k,0,1,0,1,k,0,0,0,1};
}
Matrix getB(int k){return Matrix{1,k,0,1,0,0,1,0,1,0,k,0,0,1};
}
vector<int> a,b;
struct SegmentTree{struct Node{Vec v;Matrix tag=I;bool lazy=0;};vector<Node> t;void init(int n){t=vector<Node>(n<<2);build(1,1,n);}void pushup(int p){t[p].v=t[p<<1].v+t[p<<1|1].v;}void pushtag(int p,const Matrix &m){t[p].v=t[p].v*m;t[p].tag=t[p].tag*m;t[p].lazy=1;}void pushdown(int p){if(t[p].lazy){pushtag(p<<1,t[p].tag);pushtag(p<<1|1,t[p].tag);t[p].lazy=0;t[p].tag=I;}}void build(int p,int l,int r){if(l==r){t[p].v={a[l],b[l],a[l]*b[l],0,1};return;}int mid=l+r>>1;build(p<<1,l,mid);build(p<<1|1,mid+1,r);pushup(p);}void modify(int p,int l,int r,int ql,int qr,const Matrix &x){if(ql<=l&&r<=qr){pushtag(p,x);return;}pushdown(p);int mid=l+r>>1;if(ql<=mid)modify(p<<1,l,mid,ql,qr,x);if(qr>mid)modify(p<<1|1,mid+1,r,ql,qr,x);pushup(p);}int query(int p,int l,int r,int ql,int qr){if(ql<=l&&r<=qr)return t[p].v.x03;pushdown(p);int res=0,mid=l+r>>1;if(ql<=mid)res+=query(p<<1,l,mid,ql,qr);if(qr>mid)res+=query(p<<1|1,mid+1,r,ql,qr);return res;}
} t;
void solve(){int n,m;cin>>m>>n;vector<int> sta(n+1),stb(n+1);a.assign(n+1,0);b.assign(n+1,0);for(int i=1;i<=n;i++)cin>>a[i];for(int i=1;i<=n;i++)cin>>b[i];t.init(n);cin>>m;vector<PII> Que[n+1];for(int i=0;i<m;i++){int l,r;cin>>l>>r;Que[r].emplace_back(l,i);}vector<int> ans(m);int at=0,bt=0;for(int i=1;i<=n;i++){while(at&&a[sta[at]]<a[i]){t.modify(1,1,n,sta[at-1]+1,sta[at],getA(a[i]-a[sta[at]]));--at;}sta[++at]=i;while(bt&&b[stb[bt]]<b[i]){t.modify(1,1,n,stb[bt-1]+1,stb[bt],getB(b[i]-b[stb[bt]]));--bt;}stb[++bt]=i;t.modify(1,1,n,1,i,C);for(auto [l,id]:Que[i]){ans[id]=t.query(1,1,n,l,i);}}for(int i=0;i<m;i++)cout<<ans[i]<<"\n";
}

这样就可以通过此题了，复杂度大概能快个10倍

CF1824D

首先一样的套路，离线从左到右扫，差分一下,即求${\sum }_{i=l}^r{\sum }_{j=1}^{y}g(i,j)-{\sum }_{i=l}^{x-1}g(i,j)$

用一个set维护最后一个数出现的位置，那么$g(i,j)就是set里面第一个大于等于$i$的位置；$对于一个$g(i,j)$，考虑哪些$i$的贡献发生改变，手玩一下发现是一段区间覆盖，可以用区间加代替

那么区间加，区间历史和，就是套版子了
维护向量$[sum,his,len]$表示和，历史和，长度

区间加矩阵：
$$\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ k& 0& 1\end{array}\right)$$
更新矩阵
$$\left(\begin{array}{ccc}1& 1& 0\\ 0& 1& 0\\ 1& 0& 1\end{array}\right)$$

打表简化常数，最终发现只需要维护6个变量
按照之前的步骤写，不需要手动卡常也能过

struct Matrix{int x00,x01,x11,x20,x21,x22;
};
struct Vec{int x00,x01,x02;
};
Matrix operator*(const Matrix &a,const Matrix &b){Matrix res;res.x00=a.x00*b.x00;res.x01=a.x00*b.x01+a.x01*b.x11;res.x11=a.x11*b.x11;res.x20=a.x20*b.x00+a.x22*b.x20;res.x21=a.x20*b.x01+a.x21*b.x11+a.x22*b.x21;res.x22=a.x22*b.x22;return res;
}
Vec operator * (const Vec &a,const Matrix &b){Vec res;res.x00=a.x00*b.x00+a.x02*b.x20;res.x01=a.x00*b.x01+a.x01*b.x11+a.x02*b.x21;res.x02=a.x02*b.x22;return res;
}
Vec operator +(const Vec &a,const Vec &b){return {a.x00+b.x00,a.x01+b.x01,a.x02+b.x02,};
}
Matrix I={1,0,1,0,0,1};
Matrix C={1,1,1,0,0,1};
Matrix getA(int k){return Matrix{1,0,1,k,0,1};
}
vector<int> a,b;
struct SegmentTree{struct Node{Vec v;Matrix tag=I;bool lazy=0;};vector<Node> t;void init(int n){t=vector<Node>(n<<2);build(1,1,n);}void pushup(int p){t[p].v=t[p<<1].v+t[p<<1|1].v;}void pushtag(int p,const Matrix &m){t[p].v=t[p].v*m;t[p].tag=t[p].tag*m;t[p].lazy=1;}void pushdown(int p){if(t[p].lazy){pushtag(p<<1,t[p].tag);pushtag(p<<1|1,t[p].tag);t[p].lazy=0;t[p].tag=I;}}void build(int p,int l,int r){if(l==r){t[p].v={0,0,1};return;}int mid=l+r>>1;build(p<<1,l,mid);build(p<<1|1,mid+1,r);pushup(p);}void modify(int p,int l,int r,int ql,int qr,const Matrix &x){if(ql<=l&&r<=qr){pushtag(p,x);return;}pushdown(p);int mid=l+r>>1;if(ql<=mid)modify(p<<1,l,mid,ql,qr,x);if(qr>mid)modify(p<<1|1,mid+1,r,ql,qr,x);pushup(p);}int query(int p,int l,int r,int ql,int qr){if(ql<=l&&r<=qr)return t[p].v.x01;pushdown(p);int res=0,mid=l+r>>1;if(ql<=mid)res+=query(p<<1,l,mid,ql,qr);if(qr>mid)res+=query(p<<1|1,mid+1,r,ql,qr);return res;}
} t;
void solve(){int n,q;cin>>n>>q;vector<int> a(n+1);for(int i=1;i<=n;i++)cin>>a[i];vector<int> ans(q);vector<array<int,4>> Que[n+1];for(int i=0;i<q;i++){int l,r,L,R;cin>>l>>r>>L>>R;Que[R].push_back({l,r,i,1});Que[L-1].push_back({l,r,i,-1});}t.init(n);set<int> s={0};vector<int> pre(n+1);for(int i=1;i<=n;i++){if(pre[a[i]]){auto j=s.lower_bound(pre[a[i]]);int l=*prev(j)+1,r=*j;int v=(next(j)==s.end())?i-*j:*next(j)-*j;// cout<<l<<" "<<r<<" "<<v<<"\n";t.modify(1,1,n,l,r,getA(v));s.erase(j);}t.modify(1,1,n,i,i,getA(i));pre[a[i]]=i;s.insert(i);t.modify(1,1,n,1,i,C);for(auto [l,r,id,op]:Que[i]){r=min(r,i);ans[id]+=op*t.query(1,1,n,l,r);}}for(int i=0;i<q;i++)cout<<ans[i]<<"\n";
}

P9990/2020 ICPC EcFinal G

#include<bits/stdc++.h>
#define N 1000005
#define rd read()
#define int long long
using namespace std;
int n,m,p[N],a[N],ans[N],t,l,r,op,las[N];
vector<pair<int,int> >q[N]; 
stack<int>s;
inline int read()
{int x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}while(ch>='0' && ch<='9')x=x*10+ch-'0',ch=getchar();return x*f;
}
struct ver{int x11,x21,x31;ver operator + (const ver &o){return {x11+o.x11,x21+o.x21,x31+o.x31};}
}tr[N<<2];
struct mat{int x11,x12;int x21,x22;int x31,x32;mat operator * (const mat &o){mat res;res.x11=x11*o.x11+x12*o.x21;res.x12=x11*o.x12+x12*o.x22;res.x21=x21*o.x11+x22*o.x21;res.x22=x21*o.x12+x22*o.x22;res.x31=x31*o.x11+x32*o.x21+o.x31;res.x32=x31*o.x12+x32*o.x22+o.x32;return res;}ver operator * (const ver &o){return {x11*o.x11+x12*o.x21,x21*o.x11+x22*o.x21,x31*o.x11+x32*o.x21+o.x31};}
}tag[N<<2];
inline void pushup(int k){tr[k]=tr[k<<1]+tr[k<<1|1];}
inline void build(int k,int l,int r){tag[k]={1,0,0,1,0,0};if(l==r){return tr[k].x21=1,void();}int mid=l+r>>1;build(k<<1,l,mid);build(k<<1|1,mid+1,r);pushup(k);
}
inline void add(int k,mat v){tag[k]=v*tag[k];tr[k]=v*tr[k];
}
inline void pushdown(int k){add(k<<1,tag[k]);add(k<<1|1,tag[k]);tag[k]={1,0,0,1,0,0};
}
inline void modify(int k,int l,int r,int x,int y,mat v){if(x<=l&&r<=y){return add(k,v),void();}int mid=l+r>>1;pushdown(k);if(x<=mid) modify(k<<1,l,mid,x,y,v);if(y>mid) modify(k<<1|1,mid+1,r,x,y,v);pushup(k);
}
inline ver query(int k,int l,int r,int x,int y){if(x<=l&&r<=y){return tr[k];}int mid=l+r>>1;ver res={0,0,0};pushdown(k);if(x<=mid) res=res+query(k<<1,l,mid,x,y);if(y>mid) res=res+query(k<<1|1,mid+1,r,x,y);return res;
}
signed main(){n=rd;build(1,1,n);for(int i=1;i<=n;i++) a[i]=rd;m=rd;for(int i=1;i<=m;i++) l=rd,r=rd,q[r].push_back({i,l});for(int i=1;i<=n;i++){modify(1,1,n,las[a[i]]+1,i,{0,1,1,0,0,0});modify(1,1,n,1,n,{1,0,0,1,1,0});las[a[i]]=i;for(auto [id,l]:q[i]) ans[id]=query(1,1,n,l,i).x31;}for(int i=1;i<=m;i++) cout<<ans[i]<<'\n';return 0;
}	

拿CF ECFinal的数据测试

从上到下分别是改良后的矩阵乘法(即本文章使用的)，多重标记，矩阵乘法未减少状态，以及两个普通矩阵乘法

可以看到优化还是很明显的，比赛的时候时间充足找不到打多重标记的做法可以使用，优点就是非常套路，缺点就是代码量会大一点，并且非常难调(几个初始矩阵一定要用对)