数据结构维护矩阵乘法

水博客 ing

引入：矩阵乘法

应用

题目中统计 dp 值等信息时，统计方式可以使用矩阵表达，然后多次进行统计相当于多次乘上某个矩阵。而对应的矩阵乘法（可能式广义矩阵乘法，如 $A_{i,j}=\min_k(B_{i,k}+C_{k,j})$）满足结合律时，可以使用某些方法优化矩阵积的统计，例如快速幂或者线段树优化矩阵乘法。（注意有些时候矩阵乘法不满足交换律，需要注意矩阵乘法的顺序）

数据结构（一般会是线段树）维护矩阵乘法应用较为广泛，在很多时候均有用到（某些矩阵会进行变化，使用线段树维护区间变化就很方便；然而视要求也可以 平衡树维护矩阵乘法），实现难度不等，从 ABC246H 01? Queries 到 CF1368H2 Breadboard Capacity (hard version) 都涉及到了矩阵乘法，但是某些时候修改/统计矩阵的方式不同造成了代码难度的差异。

在对区间的矩阵整体进变化时（e.g. 区间反转），可以维护区间内操作/未操作对应的矩阵积，修改时直接 swap 对应的矩阵积即可。

当然可以把矩阵乘法维护的信息放到树上。（e.g. P8820 [CSP-S 2022] 数据传输）

例题1：P7739 [NOI2021] 密码箱

解法

首先考虑以 $a_l+\frac 1{a_{l+1}+\frac1{{_{\ddots+a_{r-1}+\frac 1{a_r}}}}}$ 形式表示的 $f(l,r)$ 能不能直接展开而不进行约分。若这样写的 $f(l+1,r)$ 不用约分，则 $\frac 1{f(l+1,r)}$ 同样不用约分，$a_l+\frac1{f(l+1,r)}$ 也不用约分；如果将 $f(r,r)$ 写成 $\frac {a_r}1$ 的形式，归纳法可证明 $f(l,r)$ 可以不用约分。

考虑 $f(l,r)$ 如何快速求出。由于我们要分别计算分子分母，考虑某个 $a$ 值对 $f$ 的分子分母各自的贡献。设 $f(l+1,r)$ 的分数形式为 $\frac xy$，则 $f(l,r)$ 的分数形式为 $\frac{a_lx+y}x$，考虑将 $f$ 的分数形式写成矩阵：将 $f(l+1,r)$ 写为 $\begin{bmatrix}x\\y\end{bmatrix}$，$f(l,r)$ 写为 $\begin{bmatrix}a_lx+y\\x\end{bmatrix}$，则 $\begin{bmatrix}a_l&1\\1&0\end{bmatrix}\times\begin{bmatrix}x\\y\end{bmatrix}=\begin{bmatrix}a_lx+y\\x\end{bmatrix}$。这样计算答案就是区间矩阵积的形式（并且是从左乘到右）。

接着考虑把 $W,E$ 操作对应的矩阵 $M_W,M_E$ 写出来。由题得 $\begin{bmatrix}a&1\\1&0\end{bmatrix}\times M_W=\begin{bmatrix}a+1&1\\1&0\end{bmatrix}$，$\begin{bmatrix}a&1\\1&0\end{bmatrix}\times \begin{bmatrix}1&1\\1&0\end{bmatrix}\times M_E=\begin{bmatrix}a+1&1\\1&0\end{bmatrix}\times\begin{bmatrix}1&1\\1&0\end{bmatrix}$，$\begin{bmatrix}b+1&1\\1&0\end{bmatrix}\times M_E=\begin{bmatrix}b&1\\1&0\end{bmatrix}\times\begin{bmatrix}1&1\\1&0\end{bmatrix}\times\begin{bmatrix}1&1\\1&0\end{bmatrix}$；可得 $M_W=\begin{bmatrix}0&1\\1&1\end{bmatrix},M_E=\begin{bmatrix}2&1\\-1&0\end{bmatrix}$。

使用文艺平衡树维护矩阵积，最后计算 $\begin{bmatrix}1&0\\1&1\end{bmatrix}\times$ 矩阵积 $\times \begin{bmatrix}1\\0\end{bmatrix}$ 即可。（$\begin{bmatrix}a&1\\1&0\end{bmatrix}\times\begin{bmatrix}1\\0\end{bmatrix}=\begin{bmatrix}a\\1\end{bmatrix}$，以此避免维护 $a$ 除最后一项外的对应矩阵积和最后一项本身）

似乎有 可以不需要维护四种矩阵对应的乘积的方法。

代码

不太长。（可能是因为写过了一遍）

点此查看代码

#include <bits/stdc++.h>
using namespace std;
const int maxn=100010;
const int md=998244353;
mt19937 Rand(time(0));
int n,q,i,l,r,x,y,z,rt,tot;
char s[maxn],c;
struct mat{
    int xx,xy,yx,yy;
    inline mat operator *(const mat &a){
        return (mat){(1LL*xx*a.xx+1LL*xy*a.yx)%md,
                     (1LL*xx*a.xy+1LL*xy*a.yy)%md,
                     (1LL*yx*a.xx+1LL*yy*a.yx)%md,
                     (1LL*yx*a.xy+1LL*yy*a.yy)%md};
    }
};
const mat M[2]={{1,0,1,1},{2,1,md-1,0}};
struct node{
    int ls,rs,siz; 
    unsigned key; 
    bool rev,flp,val; 
    mat mul,rem,flm,rfm;
}tr[maxn<<1];
#define ls(p) tr[p].ls
#define rs(p) tr[p].rs
#define siz(p) tr[p].siz
#define key(p) tr[p].key
#define rev(p) tr[p].rev
#define flp(p) tr[p].flp
#define val(p) tr[p].val
#define mul(p) tr[p].mul
#define rem(p) tr[p].rem
#define flm(p) tr[p].flm
#define rfm(p) tr[p].rfm
inline int New(bool f){
    tr[++tot]={0,0,1,Rand(),0,0,f,M[f],M[!f],M[f],M[!f]};
    return tot;
}
inline void PushFlp(int p){
    if(!p) return;
    flp(p)^=1;
    swap(ls(p),rs(p));
    swap(mul(p),flm(p));
    swap(rem(p),rfm(p));
}
inline void PushRev(int p){
    if(!p) return;
    rev(p)^=1; val(p)^=1;
    swap(mul(p),rem(p));
    swap(flm(p),rfm(p));
}
inline void Pushdown(int p){
    if(flp(p)){
        PushFlp(ls(p));
        PushFlp(rs(p));
        flp(p)=0;
    }
    if(rev(p)){
        PushRev(ls(p));
        PushRev(rs(p));
        rev(p)=0;
    }
}
inline void Pushup(int p){
    siz(p)=siz(ls(p))+siz(rs(p))+1;
    mul(p)=flm(p)=M[val(p)];
    rem(p)=rfm(p)=M[!val(p)];
    if(ls(p)){
        mul(p)=mul(ls(p))*mul(p);
        rem(p)=rem(ls(p))*rem(p);
        flm(p)=flm(p)*flm(ls(p));
        rfm(p)=rfm(p)*rfm(ls(p));
    }
    if(rs(p)){
        mul(p)=mul(p)*mul(rs(p));
        rem(p)=rem(p)*rem(rs(p));
        flm(p)=flm(rs(p))*flm(p);
        rfm(p)=rfm(rs(p))*rfm(p);
    }
}
void sSplit(int p,int s,int &x,int &y){
    if(!p){
        x=y=0;
        return;
    }
    Pushdown(p);
    if(siz(ls(p))<s){
        x=p;
        sSplit(rs(p),s-siz(ls(p))-1,rs(p),y);
    }
    else{
        y=p;
        sSplit(ls(p),s,x,ls(p));
    }
    Pushup(p);
}
int Merge(int x,int y){
    if(!(x&&y)) return x|y;
    if(key(x)<key(y)){
        Pushdown(x);
        rs(x)=Merge(rs(x),y);
        Pushup(x); return x;
    }
    else{
        Pushdown(y);
        ls(y)=Merge(x,ls(y));
        Pushup(y); return y;
    }
}
inline void Ans(){
    printf("%d %d\n",mul(rt).xx,
                    (mul(rt).xx+mul(rt).yx)%md);
}
int main(){
    scanf("%d%d%s",&n,&q,s+1);
    for(i=1;i<=n;++i) rt=Merge(rt,New(s[i]=='E'));
    Ans();
    while(q--){
        scanf("%s",s+1);
        if(s[1]=='A'){
            scanf(" %c",&c);
            rt=Merge(rt,New(c=='E'));
        }
        else{
            scanf("%d%d",&l,&r);
            sSplit(rt,r,y,z);
            sSplit(y,l-1,x,y);
            if(s[1]=='F') PushRev(y);
            else PushFlp(y);
            rt=Merge(Merge(x,y),z);
        }
        Ans();
    }
    return 0;
}

例题2：P7453 [THUSCH2017] 大魔法师

解法

区间修改可看成区间乘上某个矩阵，写出矩阵乘法的形式并不难，关键在如何维护一个类似于矩阵和的形式（题目中的区间 $\sum A,\sum B,\sum C$）。

由于矩阵乘法满足分配律，则可以直接在区间上维护矩阵和，可以使用线段树维护矩阵和，还有用维护懒标记的方式维护乘上的矩阵积。