【笔记】矩阵快速幂

矩阵快速幂

矩阵乘法 + 快速幂

矩阵加法：

定义矩阵 \(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;
}