# 「学习笔记」矩阵乘法与矩阵快速幂

2022-08-05 17:34:37  阅读：11  来源： 互联网

# 「学习笔记」矩阵乘法与矩阵快速幂

## 矩阵乘

### 算法

$C_{i,j}=\sum_{k=1}^{m}A_{i,k}*B_{k,j}$

### 代码

const ll N=110,inf=1ll<<40;
ll n,m,p;
class Matrix{
public:
ll mat[N][N];
inline ll* operator[](ll x){return mat[x];}
inline Matrix operator*(Matrix ma)const{
Matrix mt;
_for(i,1,n){
_for(j,1,p){
mt[i][j]=0;
_for(k,1,m)mt[i][j]+=mat[i][k]*ma[k][j];
}
}
return mt;
}
}a,b;
inline void print(Matrix ma){
_for(i,1,n){
_for(j,1,p)printf("%lld ",ma[i][j]);
puts("");
}
return;
}
namespace SOLVE{
inline ll rnt(){
ll x=0,w=1;char c=getchar();
while(!isdigit(c)){if(c=='-')w=-1;c=getchar();}
while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x*w;
}
inline void In(){
n=rnt(),m=rnt();
_for(i,1,n)_for(j,1,m)a[i][j]=rnt();
p=rnt();
_for(i,1,m)_for(j,1,p)b[i][j]=rnt();
print(a*b);
return;
}
}


## 矩阵快速幂

### 代码（模板题）

const ll N=110,P=1e9+7,inf=1ll<<40;
ll n,k;
class Mat{
public:
ll a[N][N];
inline ll* operator[](ll x){return a[x];}
inline void one(){_for(i,1,n)a[i][i]=1;}
inline Mat operator*(Mat ma){
Mat mt;
_for(i,1,n){
_for(j,1,n){
mt[i][j]=0;
_for(k,1,n)mt[i][j]=(mt[i][j]+a[i][k]*ma[k][j]%P)%P;
}
}
return mt;
}
}a;
inline void printf(Mat ma){
_for(i,1,n){
_for(j,1,n)printf("%lld ",ma[i][j]);
puts("");
}
return;
}
inline Mat fpow(Mat ma,ll b){
Mat ans;ans.one();
while(b){
if(b&1)ans=ans*ma;
ma=ma*ma,b>>=1;
}
return ans;
}
namespace SOLVE{
inline ll rnt(){
ll x=0,w=1;char c=getchar();
while(!isdigit(c)){if(c=='-')w=-1;c=getchar();}
while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x*w;
}
inline void In(){
n=rnt(),k=rnt();
_for(i,1,n)_for(j,1,n)a[i][j]=rnt();
printf(fpow(a,k));
return;
}
}


## 练习题

### 斐波那契数列

#### 思路

\begin{aligned} F(n+1)&=\left[\begin{matrix}Fib_n+1&Fib_n\end{matrix}\right]\\ &=\left[\begin{matrix}Fib_n*1+Fib_{n-1}*1&Fib_n*1+Fib_{n-1}*0\end{matrix}\right]\\ &=\left[\begin{matrix}F(n)_{1,1}*1+F(n)_{1,2}*1&F(n)_{1,1}*1+F(n)_{1,2}*0\end{matrix}\right]\\ \end{aligned}

\begin{aligned} F(n+1)&=\left[\begin{matrix}F(n)_{1,1}\times M_{1,1}+F(n)_{1,2}\times M_{2,1}&F(n)_{1,1}\times M_{1,2}+F(n)_{1,2}\times M_{2,2}\end{matrix}\right]\\ &=F(n)\times M \end{aligned}

#### 代码

const ll N=5,inf=1ll<<40;
ll n,m;
class Mat{
public:
ll a[N][N];
inline ll* operator[](ll x){return a[x];}
friend Mat Mul(Mat m1,Mat m2){
Mat an;
_for(i,1,2){
_for(j,1,2){
an[i][j]=0;
_for(k,1,2)an[i][j]=(an[i][j]+m1[i][k]*m2[k][j]%m)%m;
}
}
return an;
}
};
inline Mat fpow(ll b){
Mat ans;ans[1][1]=ans[1][2]=1,ans[2][1]=ans[2][2]=0;
Mat ma;ma[1][1]=ma[1][2]=ma[2][1]=1,ma[2][2]=0;
while(b>0){
if(b&1)ans=Mul(ans,ma);
ma=Mul(ma,ma),b>>=1;
}
return ans;
}
namespace SOLVE{
inline ll rnt(){
ll x=0,w=1;char c=getchar();
while(!isdigit(c)){if(c=='-')w=-1;c=getchar();}
while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x*w;
}
inline void In(){
n=rnt(),m=rnt();
Mat ans=fpow(n-2);
printf("%lld\n",ans[1][1]);
return;
}
}


### [SCOI2009] 迷路

#### 思路

$f_{i,j}=\sum_{(i,k)\in\mathbb{E}}f_{k,j-w_{i,k}}$

$f_{i,j}=\sum_{1\le j\le n}w_{i,k}\times f_{k,j-1}$

#### 代码

const ll N=110,P=2009,inf=1ll<<40;
ll n,m,t,g[N][N];
class Mat{
public:
ll a[N][N];
inline ll* operator[](ll x){return a[x];}
inline void one(){_for(i,1,m)a[i][i]=1;return;}
inline Mat operator*(Mat ma){
Mat ans;
_for(i,1,m){
_for(j,1,m){
ans[i][j]=0;
_for(k,1,m)ans[i][j]=(ans[i][j]+a[i][k]*ma[k][j]%P)%P;
}
}
return ans;
}
}tu;
namespace SOLVE{
inline ll rnt(){
ll x=0,w=1;char c=getchar();
while(!isdigit(c)){if(c=='-')w=-1;c=getchar();}
while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x*w;
}
inline ll rsnt(){
char c=getchar();
while(!isdigit(c))c=getchar();
return (c^48);
}
inline Mat fpow(Mat ma,ll b){
Mat ans;ans.one();
while(b){
if(b&1)ans=ans*ma;
ma=ma*ma,b>>=1;
}
return ans;
}
inline void Zhuan(){
_for(i,1,n){
_for(j,1,9){
ll k=10*(i-1)+j;
tu[k][k+1]=1;
}
_for(j,1,n){
ll k1=10*(i-1);
ll k2=10*(j-1);
if(g[i][j])tu[k1+g[i][j]][k2+1]=1;
}
}
return;
}
inline void In(){
n=rnt(),t=rnt();
_for(i,1,n)_for(j,1,n)g[i][j]=rsnt();
m=10*n,Zhuan();
Mat ans=fpow(tu,t);
printf("%lld\n",ans[1][m-9]);
return;
}
}


### 佳佳的 Fibonacci

#### 思路

\begin{aligned} T(n)&=(F_1+F_2+\cdots+F_n)+(F_2+F_3\cdots+F_n)+(F_3+F_4+\cdots+F_n)+\cdots+(F_{n-1}+F_n)+F_n\\ &=(S(n)-S(0))+(S(n)-S(1))+(S(n)-S(3))+\cdots+(S(n)-S(n-2))+(S(n)-S(n-1))\\ &=\sum_{i=1}^{n}S(n)-S(i-1) \end{aligned}

\begin{aligned} S(n)&=\sum_{i=1}^{n}F_i\\ &=\sum_{i=1}^{n}F_{i-1}+F_{i-2}\\ &=\sum_{i=0}^{n-1}F_{i}+\sum_{i=-1}^{n-2}F_{i}\\ &=\sum_{i=0}^{n-1}F_{i}+\sum_{i=0}^{n-2}F_{i}+F_{-1}\\ &=S(n-1)+S(n-2)+1 \end{aligned}

\begin{aligned} \because &F_n=F_{n-1}+F_{n-2}\\ \therefore &F_1=F_0+F_{-1}\\ &1=0+F_{-1}\\ \therefore &F_{-1}=1\\ &同理可得：\\ &F_{-2}=-1,F_{-3}=2,F_{-4}=-3,F_{-5}=5,\cdots,F_{-n}=-1^{n+1}F_n \end{aligned}

\begin{aligned} T(n)&=\sum_{i=1}^{n}S(n)-S(i-1)\\ &=\sum_{i=1}^{n}(F_{n+2}-1)-(F_{i+1}-1)\\ &=nF_{n+2}-\sum_{i=2}^{n+1}F_{i}\\ &=nF_{n+2}-(S(n+1)-1)\\ &=nF_{n+2}-F_{n+3}+2 \end{aligned}

#### 代码

const ll N=110,inf=1ll<<40;
ll n,m;
class Mat{
public:
ll a[N][N];
inline ll* operator[](ll x){return a[x];}
inline void one(){a[1][1]=a[1][2]=1;}
inline Mat operator*(Mat ma){
Mat ans;
_for(i,1,2){
_for(j,1,2){
ans[i][j]=0;
_for(k,1,2)ans[i][j]=(ans[i][j]+a[i][k]*ma[k][j]%m)%m;
}
}
return ans;
}
};
namespace SOLVE{
inline ll rnt(){
ll x=0,w=1;char c=getchar();
while(!isdigit(c)){if(c=='-')w=-1;c=getchar();}
while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x*w;
}
inline Mat fpow(Mat ma,ll b){
Mat ans;ans.one();
while(b){
if(b&1)ans=ans*ma;
ma=ma*ma,b>>=1;
}
return ans;
}
inline void In(){
n=rnt(),m=rnt();
Mat mat;mat[1][1]=mat[1][2]=mat[2][1]=1;
Mat ans=fpow(mat,n+1);
printf("%lld\n",(n*ans[1][2]%m-ans[1][1]+m+2)%m);
return;
}
}


### 选拔队员（不知道教练从哪里找的）

#### 思路

\begin{aligned} f_{i,0}&=f_{i-1,0}+f_{i-1,1}\\ f_{i,1}&=f_{i-1,0} \end{aligned}

$\left[\begin{matrix} f_{n,0}&f_{n,1} \end{matrix}\right] \times \left[\begin{matrix} 1&1\\ 1&0 \end{matrix}\right] = \left[\begin{matrix} f_{n+1,0}&f_{n+1,1} \end{matrix}\right]$

#### 代码

const ll N=110,inf=1ll<<40;
ll T,m,n;
class Mat{
public:
ll a[5][5];
inline ll* operator[](ll x){return a[x];}
inline Mat operator*(Mat ma){
Mat ans;
_for(i,1,2){
_for(j,1,2){
ans[i][j]=0;
_for(k,1,2)ans[i][j]=(ans[i][j]+a[i][k]*ma[k][j]%m)%m;
}
}
return ans;
}
};
namespace SOLVE{
inline ll rnt(){
ll x=0,w=1;char c=getchar();
while(!isdigit(c)){if(c=='-')w=-1;c=getchar();}
while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x*w;
}
inline Mat fpow(ll b){
Mat ans;ans[1][1]=ans[1][2]=1,ans[2][1]=ans[2][2]=0;
Mat ma;ma[1][1]=ma[1][2]=ma[2][1]=1,ma[2][2]=0;
while(b){
if(b&1)ans=ans*ma;
ma=ma*ma,b>>=1;
}
return ans;
}
inline void In(){
n=rnt();
printf("%lld\n",fpow(n)[1][1]%m);
return;
}
}


### Tr A

#### 代码

const ll N=20,P=9973,inf=1ll<<40;
ll T,n,k;
class Mat{
public:
ll a[N][N];
inline ll* operator[](ll x){return a[x];}
inline void one(){_for(i,1,n)a[i][i]=1;return;}
inline void zero(){memset(a,0,sizeof(a));return;}
inline Mat operator*(Mat ma){
Mat ans;ans.zero();
_for(i,1,n){
_for(j,1,n){
ans[i][j]=0;
_for(k,1,n)ans[i][j]=(ans[i][j]+a[i][k]*ma[k][j]%P)%P;
}
}
return ans;
}
}a;
namespace SOLVE{
inline ll rnt(){
ll x=0,w=1;char c=getchar();
while(!isdigit(c)){if(c=='-')w=-1;c=getchar();}
while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x*w;
}
inline Mat fpow(Mat ma,ll b){
Mat ans;ans.zero();ans.one();
while(b){
if(b&1)ans=ans*ma;
ma=ma*ma,b>>=1;
}
return ans;
}
inline ll GetAnswer(Mat ma){
ll num=0;
_for(i,1,n)num=(num+ma[i][i]);
return num%P;
}
inline void In(){
n=rnt(),k=rnt();
a.zero();
_for(i,1,n)_for(j,1,n)a[i][j]=rnt();
Mat ans=fpow(a,k);