矩阵快速幂入门

矩阵乘法

基础概念

矩阵乘法前置数学知识：Link
大小为 \(p\) 矩阵乘法的运算规律为：

\(C_{ij}=\Sigma_{k=1}^{p}A_{ik}\times B_{kj}\)

对于两个矩阵A,B：

a1 a2   b1 b2
a3 a4   b3 b4

做乘法得到的结果矩阵C：

a1*b1+a2*b3 a1*b2+a2*b4
a3*b1+a4*b3 a3*b2+a4*b4

代码实现

	for(int k=1;k<=n;k++){
		for(int i=1;i<=n;i++){
			for(int j=1;j<=n;j++){
				res.a[i][j]=res.a[i][j]+x.a[i][k]*y.a[k][j];
			}
		}
	}

单位矩阵

注意到，形如下矩图的矩阵有着与任意相同大小的矩阵相乘，均得到原矩阵的特性。

1 0
0 1

这样的矩阵被称为单位矩阵，在处理快速幂的时候会被用到。

矩阵快速幂

因为矩阵乘法结合律成立，即存在：
\(A(BC)=(AB)C\)
因此，对于多次的矩阵乘法操作，可以进行快速幂操作，这样就可以把原来\(O(n)\)的操作转化为\(O(log(n))\)的操作，可以应对较大的数据。
将二者的操作结合在一起，用结构体来写会比较方便。

代码实现

LL n,k;
struct matrix{
	LL a[MAXN][MAXN];
	matrix(){
		memset(a,0,sizeof a);//初始化默认全为0
	}
	inline void build(){
		for(int i=1;i<=n;i++)a[i][i]=1;//单位矩阵
	}
};
matrix operator * (const matrix &x,const matrix &y){
	matrix res;
	for(int k=1;k<=n;k++){
		for(int i=1;i<=n;i++){
			for(int j=1;j<=n;j++){
				res.a[i][j]+=x.a[i][k]*y.a[k][j];
			}
		}
	}
	return res;
}//重载乘法运算符
matrix operator ^ (const matrix &x,LL &y){
	matrix res,base=x;
	res.build();
	while(y){
		if(y%2){
			res=res*base;
		}
		y/=2;
		base=base*base;
	}
	return res;
}//重载幂运算符

例题

矩阵快速幂模板普及/提高-

Link: 洛谷 P3390

就真的特么只是模板而已。

斐波那契数列普及/提高-

Link: 洛谷 P1962

根据题意，可以得到 \(f_{i}=f_{i-1}+f_{i-2}\) 的递推式。在遇到这类递推式的时候，构建矩阵是一种常见的思路。
可以构建如下的矩阵：

f[i] f[i-1]

系数矩阵：

1 0
1 1

那么该矩阵与系数矩阵相乘后，将会得到

f[i] f[i]+f[i-1]

那么达到了求斐波那契数列的效果，此时对系数矩阵先进行快速幂，再与数列的初始2项 \({1,1}\) 相乘，就可以求得答案。

代码实现

#include<bits/stdc++.h>
#define MAXN 105
#define LL long long
#define inf 0x3f3f3f3f
using namespace std;
inline LL read(){
	LL res=0,fl=1;
	char ch=getchar();
	while(!(ch>='0' && ch<='9')){if(ch=='-')fl=-1;ch=getchar();}
	while(ch>='0' && ch<='9')res=res*10+ch-'0',ch=getchar();
	return res*fl;
}
inline LL max(LL a,LL b){return a>b?a:b;}
inline LL min(LL a,LL b){return a<b?a:b;}
LL n,MOD=1e9+7,si=2;//si指的是矩阵的大小，为避免与多数题目中n混淆而换用变量名
struct matrix{
	LL a[MAXN][MAXN];
	matrix(){
		memset(a,0,sizeof a);
	}
	void build(){
		for(int i=1;i<=si;i++)a[i][i]=1;
	}
};
matrix operator *(const matrix &x,const matrix &y){
	matrix res;
	for(int k=1;k<=si;k++){
		for(int i=1;i<=si;i++){
			for(int j=1;j<=si;j++){
				res.a[i][j]=(res.a[i][j]+(x.a[i][k]*y.a[k][j])%MOD)%MOD;
			}
		}
	}
	return res;
}
matrix operator ^ (const matrix &x,LL y){
	matrix res,base=x;
	res.build();
	while(y){
		if(y%2){
			res=res*base;
		}
		y/=2;
		base=base*base;
	}
	return res;
}
int main() {
	n=read();
	if(n<=2){
		cout<<1<<'\n';
		return 0;
	}//特判数列的前2项
	matrix x,p,t;
	x.a[1][1]=1,x.a[1][2]=1;//数列的前2项
	p.a[1][1]=0,p.a[1][2]=1;
	p.a[2][1]=1,p.a[2][2]=1;
	LL tt=n-2;//tt指的是进行迭代的次数，因此需要减2
	t=p^tt;
	x=x*t;
	cout<<x.a[1][2]<<'\n';
	return 0;
}

Chino的数列普及/提高-

Link: 洛谷 P5550

总体思路与上一题是类似的。在已经得到固定递推式的时候，只需要构建矩阵与系数矩阵，先进行一组操作，再对已经操作一组的矩阵进行快速幂，就可以得到 \(k\) 次操作后的结果。

对于交换 \(x,y\) 的操作，可以构建 \(a[x][y]=a[y][x]=1\) 的系数矩阵。
例如：交换 \(2,5\) 位置上的数，即构建如下系数矩阵：

同理，对于将数向前移动一位的操作，需要构建 \(a[i][i-1]=1\) 的矩阵。

代码实现

#include<bits/stdc++.h>
#define MAXN 105
#define LL long long
#define inf 0x3f3f3f3f
using namespace std;
inline LL read(){
	LL res=0,fl=1;
	char ch=getchar();
	while(!(ch>='0' && ch<='9')){if(ch=='-')fl=-1;ch=getchar();}
	while(ch>='0' && ch<='9')res=res*10+ch-'0',ch=getchar();
	return res*fl;
}
inline LL max(LL a,LL b){return a>b?a:b;}
inline LL min(LL a,LL b){return a<b?a:b;}
LL n,s,m,k,si=85;
struct matrix{
	LL a[MAXN][MAXN];
	matrix(){
		memset(a,0,sizeof a);
	}
	void build(){
		for(int i=1;i<=si;i++)a[i][i]=1;
	}
};
matrix operator *(const matrix &x,const matrix &y){
	matrix res;
	for(int k=1;k<=si;k++){
		for(int i=1;i<=si;i++){
			for(int j=1;j<=si;j++){
				res.a[i][j]+=x.a[i][k]*y.a[k][j];
			}
		}
	}
	return res;
}
matrix operator ^ (const matrix &x,LL y){
	matrix res,base=x;
	res.build();
	while(y){
		if(y%2){
			res=res*base;
		}
		y/=2;
		base=base*base;
	}
	return res;
}
int main() {
	n=read(),s=read(),m=read(),k=read();
	si=n;
	matrix opt1,opt2,ans,p;
	
	for(int i=1;i<=n;i++)ans.a[1][i]=read();
	opt1.build();
	opt1.a[s][s]=opt1.a[m][m]=0;
	opt1.a[s][m]=opt1.a[m][s]=1;
	
	for(int i=2;i<=n;i++)opt2.a[i][i-1]=1;
	opt2.a[1][n]=1;
	
	p=opt1*opt2;
	p=p^k;
	ans=ans*p;
	
	for(int i=1;i<=n;i++)cout<<ans.a[1][i]<<' ';
	return 0;
}

可乐普及+/提高

Link: 洛谷 P3758、洛谷 P5789（数据加强版）

这道题就涉及到了领接矩阵的幂的含义。
设现在存在一个领接矩阵， \(x,y\) 间的连边用 \(a[x][y]=a[y][x]=1\) 表示，那么根据 \(Floyd\) 考虑，\(A^k\) 的含义即为在 \((i,j)\) 表示从 \(i\) 到 \(j\) 走过 \(k\) 步的方案数。

有了这样的概念后，普通的连边问题就解决了。接下来考虑原地走与自爆的情况。
原地走很好考虑，直接建一条自己连自己的环就行。即 \(a[i][i]=1\) 。
自爆的情况下，可以考虑为走到一个节点后，再无法去到其他节点，那么可以建立一个 \(0\) 节点，只建立其他点到它的单向边。即 \(a[i][0]=1\) 。

代码实现

#include<bits/stdc++.h>
#define MAXN 110
#define LL long long
#define inf 0x3f3f3f3f
#define MOD 2017
using namespace std;
inline LL read(){
	LL res=0,fl=1;
	char ch=getchar();
	while(!(ch>='0' && ch<='9')){if(ch=='-')fl=-1;ch=getchar();}
	while(ch>='0' && ch<='9')res=res*10+ch-'0',ch=getchar();
	return res*fl;
}
inline LL max(LL a,LL b){return a>b?a:b;}
inline LL min(LL a,LL b){return a<b?a:b;}
LL n,m,t,si=105;
struct matrix{
	LL a[MAXN][MAXN];
	matrix(){
		memset(a,0,sizeof a);
	}
	void build(){
		for(int i=0;i<=si;i++)a[i][i]=1;
	}
};
matrix operator *(const matrix &x,const matrix &y){
	matrix res;
	for(int k=0;k<=si;k++){
		for(int i=0;i<=si;i++){
			for(int j=0;j<=si;j++){
				res.a[i][j]=(res.a[i][j]+(x.a[i][k]*y.a[k][j])%MOD)%MOD;
			}
		}
	}
	return res;
}
matrix operator ^ (const matrix &x,LL y){
	matrix res,base=x;
	res.build();
	while(y){
		if(y%2){
			res=res*base;
		}
		y/=2;
		base=base*base;
	}
	return res;
}
int main() {
	LL a,b,sum=0;
	n=read(),m=read();
	matrix x,ans;
	for(int i=1;i<=m;i++){
		a=read(),b=read();
		x.a[a][b]=1,x.a[b][a]=1;
	}
	t=read();
	for(int i=0;i<=n;i++)x.a[i][i]=1;
	for(int i=1;i<=n;i++)x.a[i][0]=1;
	ans=x^t;
	for(int i=0;i<=n;i++)sum=(sum+ans.a[1][i])%MOD;
	cout<<sum<<'\n';
	return 0;
}

posted @ 2022-06-05 11:31 Powerless233 阅读(165) 评论(0) 收藏举报

刷新页面返回顶部

Powerless's Blog

矩阵快速幂入门

矩阵乘法

基础概念

代码实现

单位矩阵

矩阵快速幂

代码实现

例题

矩阵快速幂模板 普及/提高-

斐波那契数列 普及/提高-

代码实现

Chino的数列 普及/提高-

代码实现

可乐 普及+/提高

代码实现

矩阵快速幂模板普及/提高-

斐波那契数列普及/提高-

Chino的数列普及/提高-

可乐普及+/提高