矩阵乘法与快速幂优化-优快云博客

本文链接：https://blog.youkuaiyun.com/YYHS_WSF/article/details/84930896

https://vjudge.net/contest/274056#overview

矩阵乘法：利用矩阵优化其次的递推式

矩阵和矩阵乘法相必大家都知道（不会的右转百度）

把 $n*m$ 与 $m*p$ 的矩阵相乘，得到一个 $n*p$ 的矩阵即是矩阵乘法，

时间最暴躁的为 $O(n^{3})$

———————————————————————————————————————————————————————

#include<cstdio>
using namespace std;
const int N=105;
int n,m,p,a[N][N],b[N][N],ans;
int main()
{
	scanf("%d%d",&n,&m);
	for(int i=1;i<=n;i++)
		for(int j=1;j<=m;j++)
			scanf("%d",&a[i][j]);
	scanf("%d",&p);
	for(int i=1;i<=m;i++)
		for(int j=1;j<=p;j++)
			scanf("%d",&b[i][j]);
	for(int i=1;i<=n;i++)
		for(int j=1;j<=p;j++)
		{
			ans=0;
			for(int k=1;k<=m;k++)
				ans+=a[i][k]*b[k][j];
			printf("%d%c",ans,j==p?'\n':' ');
		}
	return 0;
}

————————————————————————————————————————————————————————

矩阵乘法满足结合律，但不满足交换律（证明略）

由于它满足结合律，所以可以用快速幂优化（前提一定要是方阵（ $n*n$ 的矩阵）O_O）

记 $E$ 为单位矩阵，就是主对角线都为1，其余为0的方阵

————————————————————————————————————————————————————————

例：Fibonacci 第 n 项

递推式： $f_{i}=f_{i-1}+f_{i-2}$

$\begin{bmatrix} f_{i}\\f_{i-1} \end{bmatrix} = \begin{bmatrix} 1&1\\ 1&0 \end{bmatrix} \begin{bmatrix} f_{i-1}\\f_{i-2} \end{bmatrix}$

$\begin{bmatrix} f_{i}\\f_{i-1} \end{bmatrix} = \begin{bmatrix} 1 & 1\\ 1& 0 \end{bmatrix} ^_{i-1}$ $\begin{bmatrix} f_{1}\\f_{0} \end{bmatrix}$

矩阵快速幂即可

#include<cstdio>
#define ll long long
using namespace std;

int n,p;
struct A { ll x[10][10]; }a,ret;

inline A time(A x,A y)
{
	A z=x; 
	for(int i=1;i<=2;i++)
		for(int j=1;j<=2;j++) 
		{
			x.x[i][j]=0;
			for(int k=1;k<=2;k++)
				x.x[i][j]=(z.x[i][k]*y.x[k][j]%p+x.x[i][j])%p;
		}
	return x; 
}

void ksm(int b)
{
	ret.x[1][1]=1,ret.x[1][2]=0,
	ret.x[2][1]=0,ret.x[2][2]=1; //E(单位矩阵) 
	while(b)
	{
		if(b&1) ret=time(ret,a);
		a=time(a,a); b>>=1;
	}
}

int main()
{
	scanf("%d%d",&n,&p);
	a.x[1][1]=1,a.x[1][2]=1,
	a.x[2][1]=1,a.x[2][2]=0;
	ksm(n-1);
	printf("%d\n",ret.x[1][1]);
	return 0;
}

练：Fibonacci

当n=0是要特判

例：Fibonacci 第 n 项和

再加一维 $S$

#include<cstdio>
#define ll long long
using namespace std;

int n,p;
struct A { ll x[10][10]; }a,ret;

inline A time(A x,A y)
{
	A z=x; 
	for(int i=1;i<=3;i++)
		for(int j=1;j<=3;j++) 
		{
			x.x[i][j]=0;
			for(int k=1;k<=3;k++)
				x.x[i][j]=(z.x[i][k]*y.x[k][j]%p+x.x[i][j])%p;
		}
	return x;
}

void ksm(int b)
{
	ret.x[1][1]=1,ret.x[1][2]=0,ret.x[1][3]=0;
	ret.x[2][1]=0,ret.x[2][2]=1,ret.x[2][3]=0,
	ret.x[3][1]=0,ret.x[3][2]=0,ret.x[3][3]=1;
	while(b)
	{
		if(b&1) ret=time(ret,a);
		a=time(a,a); b>>=1;
	}
}

int main()
{
	scanf("%d%d",&n,&p);
	a.x[1][1]=1,a.x[1][2]=1,a.x[1][3]=0,
	a.x[2][1]=1,a.x[2][2]=0,a.x[2][3]=0;
	a.x[3][1]=1,a.x[3][2]=1,a.x[3][3]=1;
	ksm(n-1);
	printf("%d\n",(ret.x[3][1]+ret.x[3][3])%p);
	return 0;
}

例：(F1 +2F2+3F3 +...+nFn )modm

这类问题，往往令 $g_{i}=f_{i}*i$ ，再构造矩阵 $5*5$

#include<cstdio>
#define ll long long
using namespace std;

int n,p;
struct A { ll x[10][10]; }a,ret;

inline A time(A x,A y)
{
	A z=x; 
	for(int i=1;i<=5;i++)
		for(int j=1;j<=5;j++) 
		{
			x.x[i][j]=0;
			for(int k=1;k<=5;k++)
				x.x[i][j]=(z.x[i][k]*y.x[k][j]%p+x.x[i][j])%p;
		}
	return x;
}

void ksm(int b)
{
	for(int i=1;i<=5;i++)
		for(int j=1;j<=5;j++)
			ret.x[i][j]=i==j;
	while(b)
	{
		if(b&1) ret=time(ret,a);
		a=time(a,a); b>>=1;
	}
}

int main()
{
	scanf("%d%d",&n,&p);
	a.x[1][1]=1,a.x[1][2]=1,a.x[1][3]=0,a.x[1][4]=0,a.x[1][5]=0;
	a.x[2][1]=1,a.x[2][2]=0,a.x[2][3]=0,a.x[2][4]=0,a.x[2][5]=0;
	a.x[3][1]=0,a.x[3][2]=0,a.x[3][3]=0,a.x[3][4]=1,a.x[3][5]=0;
	a.x[4][1]=1,a.x[4][2]=2,a.x[4][3]=1,a.x[4][4]=1,a.x[4][5]=0;
	a.x[5][1]=1,a.x[5][2]=2,a.x[5][3]=1,a.x[5][4]=1,a.x[5][5]=1;
	ksm(n-1);
	printf("%lld\n",(ret.x[5][1]+ret.x[5][4]+ret.x[5][5])%p);
	return 0;
}

用途：优化DP转移

例：zoj3497:Mistwald

以邻接表为矩阵，用0/1表示可不可以走

#include<cstdio>
using namespace std;

const int N=30;
int T,n,m,c,k,TT;
inline int id(int x,int y){ return (x-1)*m+y; }
struct A{int x[N][N]; }ret,g,a;

A time(A x,A y)
{
	A z=x;
	for(int i=1;i<=c;i++)
		for(int j=1;j<=c;j++)
		{
			x.x[i][j]=0;
			for(int k=1;k<=c;k++)
			{
				x.x[i][j]|=z.x[i][k]&y.x[k][j];
				if(x.x[i][j]) break;  //无聊的剪枝 
			}
		}
	return x;
}

void ksm(int b)
{
	for(int i=1;i<=c;i++)
		for(int j=1;j<=c;j++)
			ret.x[i][j]=i==j;
	while(b)
	{
		if(b&1) ret=time(ret,a);
		b>>=1; a=time(a,a);
	}
}

int read()
{
	int ret=0; char ch=getchar();
	while(ch<'0'||ch>'9') ch=getchar();
	while(ch>='0'&&ch<='9')
		ret=(ret<<1)+(ret<<3)+ch-'0',ch=getchar();
	return ret;
}

int main()
{
	T=read();
	while(T--)
	{
		n=read(),m=read(),c=n*m;
		for(int i=1;i<=n;i++)
			for(int j=1;j<=m;j++)
			{
				int x=read(),y=read();
				if(i!=n||j!=m) g.x[id(x,y)][id(i,j)]=1;
				x=read(),y=read();
				if(i!=n||j!=m) g.x[id(x,y)][id(i,j)]=1;
				x=read(),y=read();
				if(i!=n||j!=m) g.x[id(x,y)][id(i,j)]=1;
				x=read(),y=read();
				if(i!=n||j!=m) g.x[id(x,y)][id(i,j)]=1;
			} //到(n,m)结束游戏，GG，所以不用加 
		TT=read();
		while(TT--)
		{
			k=read();
			a=g;
			ksm(k);
			if(!ret.x[c][1]) puts("False");
				else 
				{
					bool f=0;
					for(int i=1;i<c;i++)
						if(ret.x[i][1]) 
						{
							f=1; break;
						}
					puts(!f?"True":"Maybe");
				}
		}
		puts("");
	}
	return 0;
}

例：Kiki & Little Kiki 2

双数次为不用变，单数次为要变，所以可以异或来搞一搞，

转移矩阵：第i行，第i-1项和第i项都为1，（i-1=0的话，则为n），其余都为0

#include<cstdio>
#include<cstring>
using namespace std;

const int N=105;
int n,k;
bool ans[N],f[N];
struct A{bool x[N][N]; }ret,a;
char s[N];

inline A time(A x,A y)
{
	A z=x;
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
		{
			x.x[i][j]=0;
			for(int k=1;k<=n;k++)
				x.x[i][j]^=z.x[i][k]&y.x[k][j];
		}
	return x;
}
void ksm(int b)
{
	for(int i=1;i<=n;i++) 
		for(int j=1;j<=n;j++)
			ret.x[i][j]=i==j;
	while(b)
	{
		if(b&1) ret=time(ret,a);
		a=time(a,a),b>>=1;
	}
}

int main()
{
	while(scanf("%d",&k)!=EOF)
	{
		scanf("%s",s+1);
		n=strlen(s+1);
		for(int i=1;i<=n;i++)
			for(int j=1;j<=n;j++)
				a.x[i][j]=i==j||i==1&&j==n||i==j+1;
		ksm(k); 
		for(int i=1;i<=n;i++)
			f[i]=s[i]=='1',
			ans[i]=0;
		for(int i=1;i<=n;i++)
			for(int j=1;j<=n;j++)
				ans[i]^=ret.x[i][j]&f[j];
		for(int i=1;i<=n;i++) 
			printf("%d",ans[i]);
		puts("");
	}
	return 0;
}