（简记）矩阵快速幂 矩阵加速递推 动态 DP
矩阵乘法和矩阵幂
一个 \(m \times n\) 的矩阵是一个由 \(m\) 行 \(n\) 列元素排列成的矩形阵列。即形如
本题中认为矩阵中的元素 \(a_{i j}\) 是整数。
两个大小分别为 \(m \times n\) 和 \(n \times p\) 的矩阵 \(A, B\) 相乘的结果为一个大小为 \(m \times p\) 的矩阵。将结果矩阵记作 \(C\)，则
而如果 \(A\) 的列数与 \(B\) 的行数不相等，则无法进行乘法。
可以验证，矩阵乘法满足结合律，即 \((A B) C = A (B C)\)。
一个大小为 \(n \times n\) 的矩阵 \(A\) 可以与自身进行乘法，得到的仍是大小为 \(n \times n\) 的矩阵，记作 \(A^2 = A \times A\)。进一步地，还可以递归地定义任意高次方 \(A^k = A \times A^{k - 1}\)，或称 \(A^k = \underbrace{A \times A \times \cdots \times A}_{k \text{ 次}}\)。
特殊地，定义 \(A^0\) 为单位矩阵 \(I = \begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{bmatrix}\)。
为什么要用到矩阵快速幂
当需要利用矩阵解决一些问题时，矩阵快速幂显得尤其好用。譬如当你知道一个序列的每个数与前若干个数的关系时，求出数列第\(k\)个数值（\(k\le 10^{18}\)）对某个数取模时。
可以帮助你在\(O(n^3\log k)\)的时间内轻松解决问题。
具体实现（代码）
只要掌握一些operator与结构体相关知识，实现并不困难，具体过程如下。
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=105,MOD=1e9+7;
int n,k;
struct matrix{
int row,col;
long long data[N][N];
matrix(int r,int c,bool isI){
row=r;
col=c;
memset(data,0,sizeof(data));
if(isI)
for(int i=0;i<row;i++)
data[i][i]=1;
}
};
matrix operator * (const matrix &a,const matrix &b){
matrix c(a.row,b.col,false);
for(int i=0;i<a.row;i++)
for(int j=0;j<b.col;j++)
for(int k=0;k<a.col;k++)
c.data[i][j]=(c.data[i][j]+a.data[i][k]*b.data[k][j])%MOD;
return c;
}
matrix powMatrix(matrix a,int k){
matrix ans(a.row,a.col,true);
for(;k;k>>=1){
if(k&1)ans=ans*a;
a=a*a;
}
return ans;
}
signed main(){
ios::sync_with_stdio(0);
cin.tie(0);cout.tie(0);
cin>>n>>k;
matrix fa(n,n,false);
for(int i=0;i<n;i++){
for(int j=0;j<n;j++){
cin>>fa.data[i][j];
}
}
fa=powMatrix(fa,k);
for(int i=0;i<n;i++){
for(int j=0;j<n;j++){
cout<<fa.data[i][j]<<' ';
}
cout<<'\n';
}
return 0;
}
构造矩阵
重定义 动态 DP
矩阵乘法的普适定义是这样的：\(C_{i,j}=\sum A_{i,k}\times B_{k,j}\)，其中应该满足 \(A\) 的列数与 \(B\) 的行数一致。
在此基础上，我们可以改变矩阵乘法原有定义，如把 \(\times\) 改成 \(+\)，使其仍然具有结合律等普遍性质。这种改变方法在P4719 【模板】动态 DP这类问题中尤其常见。在动态 DP 问题中，你需要在树上动态维护一个树形 DP，其中需要涉及到对轻重儿子的操作与区分。
对于本题矩阵方面的剖析，由于树链剖分的结果是区间左端点（链顶）的 DFS 序小，右端点（链尾）的 DFS 序大，所以需要改变以往的矩阵乘法顺序。
具体地，原矩阵长这样：
重新构造顺序后，矩阵应该长这样：
具体来说，原构造矩阵 \(base\) 中 \(base_{i,k}\) 指从 \(i\) 转移到 \(k\) 的权值。新构造矩阵 \(base\) 中的 \(base_{i,k}\) 指的是从 \(k\) 转移到 \(i\) 的权值。
运用这个转化方式，我们回到题目。
本题要求维护树形 DP 的最大独立集信息，那么可以每次 update 记录从节点 \(u\) 到链顶 \(top[u]\) 的矩阵累积，然后由于链顶一定是父亲的轻儿子，所以需要更新 \(fa[top[u]]\) 的相应信息。
矩阵加速递推
P1939 矩阵加速（数列）
题目描述
已知一个数列 \(a\)，它满足：
求 \(a\) 数列的第 \(n\) 项对 \(10^9+7\) 取余的值。
此处我们需要一个\(1\times 3\)的矩阵，形如：
由题，初始矩阵即为：
构造一个\(3\times 3\)的矩阵，形如：
其中第一列表示对第一个数的转移，第二列表示对第二个数的转移，以此类推。
具体地，如\((1,3)\)和\((3,3)\)表示进行乘法后的第\(3\)的位置受到原矩阵第\(1\)和\(3\)位置的影响，即为两者的和。
同理，如\((3,2)\)表示第\(2\)个位置为原矩阵第\(3\)为的数值。
如此便可构造出相应的矩阵。
P6435 「EZEC-1」数列
无法做出 key observation，于是无法转化出来。
需要注意到的是第 \(i\) 层的公差是 \((a+b)^{i-1}\)，然后就可以得到每行的 \(A_1,A_2\)，且下一行的 \(A_1'\) 仅跟 \(A_1,A_2\) 有关，令 \(f_i\) 表示第 \(i+1\) 行的第一项，那么依此就可以列出递推式：
矩阵加速即可。
#include<bits/stdc++.h>
using namespace std;
typedef __int128 LL;
LL n,a,b,c,MOD;
struct mat{int r,c;LL d[3][3];};
mat calc,f,base;
const mat operator *(const mat &a,const mat &b){
calc.r=a.r,calc.c=b.c;
for(int i=0;i<a.r;i++)
for(int j=0;j<b.c;j++){
calc.d[i][j]=0;
for(int k=0;k<a.c;k++)
calc.d[i][j]=(calc.d[i][j]+a.d[i][k]%MOD*(b.d[k][j]%MOD)%MOD)%MOD;
}
return calc;
}
mat qkpow(mat x,LL y){
mat res;
res.r=x.r;
res.c=x.c;
for(int i=0;i<x.r;i++)
for(int j=0;j<x.c;j++){
if(i==j)res.d[i][j]=1;
else res.d[i][j]=0;
}
while(y){
if(y&1)res=res*x;
x=x*x;
y>>=1;
}
return res;
}
int main(){
long long N,A,B,C,mod;
cin>>N>>A>>B>>C>>mod;
n=N,a=A,b=B,c=C,MOD=mod;
f.r=1,f.c=3;
base.r=3,base.c=3;
f.d[0][0]=1,f.d[0][1]=1,f.d[0][2]=1;
base.d[0][0]=a+b,base.d[1][0]=b;
base.d[1][1]=a+b,base.d[2][0]=c;
base.d[2][2]=1;
f=f*qkpow(base,n-1);
cout<<(long long)f.d[0][0];
return 0;
}
P10581 [蓝桥杯 2024 国 A] 重复的串
令 \(f_{i,e,len}\) 表示 \([1,i]\) 匹配 \(e\) 次匹配到 \(S\) 中第 \(len\) 个的方案数。加速递推 \(i\) 维，每次枚举下一个要取什么字母即可，如果取满了 \(len=|S|\)，就转移到 \(f_{i+1,e+1,nxt_{len}}\)，其中 \(nxt_i\) 数组就是 KMP 的 \(\texttt{Fail}\) 数组。无法存储两维，乘开存到一维即可。
点击查看代码
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N=105;
const LL MOD=998244353;
struct mat{
int r,c;
LL a[N][N];
mat operator *(const mat&b)const {
mat res;res.r=r;res.c=b.c;
for(int i=0;i<r;i++)
for(int j=0;j<b.c;j++)res.a[i][j]=0;
for(int i=0;i<r;i++)
for(int j=0;j<b.c;j++)
for(int k=0;k<c;k++)
(res.a[i][j]+=a[i][k]*b.a[k][j]%MOD)%=MOD;
return res;
}
}now,per;
int n,m,nxt[N];
char s[N];
inline int pack(int x,int y){return x*m+y;}
mat qkpow(mat x,int y){
mat res=x;
for(int i=0;i<res.r;i++)
for(int j=0;j<res.c;j++){
if(i==j)res.a[i][j]=1;
else res.a[i][j]=0;
}
while(y){
if(y&1)res=res*x;
x=x*x;
y>>=1;
}
return res;
}
int main(){
ios::sync_with_stdio(0);
cin.tie(0);cout.tie(0);
cin>>(s+1)>>n;m=strlen(s+1);
int pos=0;
for(int i=2;i<=m;i++){
while(pos&&s[pos+1]!=s[i])pos=nxt[pos];
if(s[pos+1]==s[i])nxt[i]=++pos;
else nxt[i]=0;
}
now.r=1;now.c=m*3;
now.a[0][0]=1;
per.r=m*3;per.c=m*3;
for(int e=0;e<=2;e++){
for(int i=0;i<m;i++){
for(char ch='a';ch<='z';ch++){
pos=i;
while(pos&&s[pos+1]!=ch)pos=nxt[pos];
pos+=(s[pos+1]==ch);
if(pos==m){
if(e==2)continue;
per.a[pack(e,i)][pack(e+1,nxt[m])]++;
}
else per.a[pack(e,i)][pack(e,pos)]++;
}
}
}
now=now*qkpow(per,n);
LL res=0;
for(int i=0;i<m;i++)
res=(res+now.a[0][pack(2,i)])%MOD;
cout<<res;
return 0;
}
一些trick
如上这类矩阵问题内，如果存在序列推导的常数项或者指数项，可以相应地加入内容。
具体内容有待补充。
浙公网安备 33010602011771号