矩阵快速幂
矩阵乘法 + 快速幂
矩阵加法：
定义矩阵 \(C=A+B\)。
\(C_{i,j}=A_{i,j}+B_{i,j}\)。
矩阵乘法：
计算两个矩阵的乘法。\(n \times m\) 阶的矩阵 \(A\) 乘以 \(m \times k\) 阶的矩阵 \(B\) 得到的矩阵 \(C\) 是 \(n \times k\) 阶的，且 \(C[i][j]=A[i][0] \times B[0][j]+A[i][1] \times B[1][j]+\) …… \(+A[i][m-1] \times B[m-1][j](C[i][j]\) 表示 \(C\) 矩阵中第 \(i\) 行第 \(j\) 列元素）。
（摘自 Luogu B2105 题面）。
即：第 \(i\) 行与第 \(j\) 列对应乘起来加和放在 \((i,j)\)。
第一个矩阵为 \(n \times m\)，第二个矩阵式 \(m \times k\)，则两个矩阵相乘为：
for(int i=1;i<=n;i++){for(int j=1;j<=k;j++){for(int d=1;d<=m;d++){c[i][j]+=a[i][d]*b[d][j];}}
}
\(O(n^3)\)。
Luogu P3390【模板】矩阵快速幂
把快速幂的乘法直接换成矩阵乘法即可。
但是要注意，结果的初始值要定义为对角线全是 \(1\) 的矩阵（单位矩阵）。单位矩阵乘任何矩阵都等于那个矩阵。
#include<iostream>
#include<cstring>
#define int long long
using namespace std;
const int p=1e9+7;
int n,k;
int a[110][110],s[110][110],c[110][110];
void sta(){memset(c,0,sizeof(c));for(int i=1;i<=n;i++){for(int j=1;j<=n;j++){for(int d=1;d<=n;d++){c[i][j]=(c[i][j]+(s[i][d]*a[d][j])%p)%p;}}}for(int i=1;i<=n;i++){for(int j=1;j<=n;j++) s[i][j]=c[i][j];}
}
void ata(){memset(c,0,sizeof(c));for(int i=1;i<=n;i++){for(int j=1;j<=n;j++){for(int d=1;d<=n;d++){c[i][j]=(c[i][j]+(a[i][d]*a[d][j])%p)%p;}}}for(int i=1;i<=n;i++){for(int j=1;j<=n;j++) a[i][j]=c[i][j];}
}
signed main(){cin>>n>>k;for(int i=1;i<=n;i++) s[i][i]=1;for(int i=1;i<=n;i++){for(int j=1;j<=n;j++){cin>>a[i][j];}}while(k){if(k&1) sta();ata();k>>=1;}for(int i=1;i<=n;i++){for(int j=1;j<=n;j++) cout<<s[i][j]<<" ";cout<<endl;}return 0;
}
运算符重载：把乘号重载成矩阵乘法的形式。
方法：把矩阵都定义成结构体，然后写运算符重载函数。
struct Matrix{int mx[110][110];
}a,s,c;
Matrix operator *(const Matrix &a,const Matrix &b){memset(c.mx,0,sizeof(c.mx));for(int i=1;i<=n;i++){for(int j=1;j<=n;j++){for(int d=1;d<=n;d++){c.mx[i][j]=(c.mx[i][j]+(a.mx[i][d]*b.mx[d][j])%p)%p;}}}return c;
}
这样快速幂就能写成：
while(k){if(k&1) s=s*a;a=a*a;k>>=1;}
同时删去了两个重复的函数及来回赋值的过程，优化了时间常数。
矩阵加速递推
Luogu P1962 斐波那契数列
可以用矩阵快速幂来加速斐波那契数列的递推。
推一个式子，
\[\begin{bmatrix} f_n & f_{n-1} \end{bmatrix} = \begin{bmatrix} f_{n-1} & f_{n-2} \end{bmatrix} \times \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}
\]
#include<iostream>
#include<cstring>
#define int long long
using namespace std;
const int p=1e9+7;
struct Matrix{int mx[3][3];
}a,s,c;
int n;
Matrix operator *(const Matrix &a,const Matrix &b){memset(c.mx,0,sizeof(c.mx));for(int i=1;i<=2;i++){for(int j=1;j<=2;j++){for(int d=1;d<=2;d++){c.mx[i][j]=(c.mx[i][j]+(a.mx[i][d]*b.mx[d][j])%p)%p;}}}return c;
}
signed main(){cin>>n;n-=1;s.mx[1][1]=s.mx[1][2]=s.mx[2][1]=1;s.mx[2][2]=0;a.mx[1][1]=a.mx[1][2]=1;while(n){if(n&1) a=s*a;s=s*s;n>>=1;}cout<<a.mx[1][1];return 0;
}
Luogu P5175 数列
#include<iostream>
#include<cstdio>
#include<cstring>
#define int long long
using namespace std;
const int p=1e9+7;
struct Matrix{int mx[5][5];
}a,s,c;
int T,n,k,f1,f2,x,y;
Matrix operator *(const Matrix &a,const Matrix &b){memset(c.mx,0,sizeof(c.mx));for(int i=1;i<=4;i++){for(int j=1;j<=4;j++){for(int d=1;d<=4;d++){c.mx[i][j]=(c.mx[i][j]%p+(a.mx[i][d]*b.mx[d][j])%p)%p;}}}return c;
}
signed main(){cin>>T;while(T--){memset(s.mx,0,sizeof(s.mx));memset(a.mx,0,sizeof(a.mx));scanf("%lld%lld%lld%lld%lld",&n,&f1,&f2,&x,&y);f1%=p;f2%=p;x%=p;y%=p;for(int i=1;i<=4;i++) s.mx[i][i]=1;a.mx[1][1]=a.mx[2][3]=1;a.mx[2][1]=a.mx[2][2]=x*x%p;a.mx[3][1]=a.mx[3][2]=y*y%p;a.mx[4][1]=a.mx[4][2]=2*x*y%p;a.mx[2][4]=x;a.mx[4][4]=y%p;if(n==1) printf("%lld\n",(f1*f1)%p);else if(n==2) printf("%lld\n",(f1*f1%p+f2*f2%p)%p);else{int s2=(f1*f1%p+f2*f2%p)%p;int f22=(f2*f2)%p;int f11=(f1*f1)%p;int f12=(f1*f2)%p;n-=2;while(n){if(n&1) s=s*a;a=a*a;n>>=1;}printf("%lld\n",(((s2*s.mx[1][1])%p+(f22*s.mx[2][1]))%p+((f11*s.mx[3][1])%p+(f12*s.mx[4][1])%p)%p)%p);}}return 0;
}
Luogu P2233 [HNOI2002] 公交车路线
构造每个站台之间的邻接矩阵，把矩阵平方，就能找到哪些站台到其他站台能换2次。把这个矩阵乘 \(n\) 次然后输出 \([1][5]\)。
注意构造矩阵时，E站台不能有出边。
运用：乘法原理/加法原理。
#include<iostream>
#include<cstring>
using namespace std;
const int dby=1000;
struct Matrix{int mx[10][10];
}c,s,a;
Matrix operator *(const Matrix &a,const Matrix &b){memset(c.mx,0,sizeof(c.mx));for(int i=1;i<=8;i++){for(int j=1;j<=8;j++){for(int d=1;d<=8;d++){c.mx[i][j]=(c.mx[i][j]%dby+(a.mx[i][d]*b.mx[d][j])%dby)%dby;}}}return c;
}
int n;
int main(){cin>>n;for(int i=1;i<=8;i++) s.mx[i][i]=1;a.mx[1][2]=a.mx[2][1]=a.mx[2][3]=a.mx[3][2]=a.mx[3][4]=a.mx[4][3]=a.mx[4][5]=a.mx[6][5]=a.mx[6][7]=a.mx[7][6]=a.mx[7][8]=a.mx[8][7]=a.mx[8][1]=a.mx[1][8]=1;while(n){if(n&1) s=s*a;a=a*a;n>>=1;}printf("%d",s.mx[1][5]%dby);return 0;
}
