前言

一般解决普通的区间历史和，只需要定义辅助$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的数据测试

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

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