#创作计划#矩阵快速幂浅析及应用

暑假神（开学祭

2024-07-03 18:24:53

发布于：上海

130阅读

0回复

0点赞

矩阵快速幂，听起来很难，模板题也都是绿色的，令人望而却步。我在发呆时突然想起了矩阵快速幂，于是花了一个下午去学习它。于是，就有了这篇文章。

ACGO 主题库上似乎还没有矩阵快速幂的题目，我从洛谷上搬运了几道。数据自造，强度应该是够的，如果搬的题目或数据有什么问题还请联系我。题单：矩阵模板题。原题分别为洛谷 P3390，P1939 和 P1962。

正文从此开始。

前置知识：

概念

先从矩阵的乘法讲起。矩阵乘法 $AB=C$ 的定义为：

$c_{i j} = \sum_{k = 1}^{n} a_{i k} b_{k j} \text{,\qquad($1 \le i \le m$, $1 \le j \le p$).}$

其中 $A$ 的大小为 $m\times n$ ， $B$ 的大小为 $n\times p$ ，结果 $C$ 的大小为 $m\times p$ 。如果 $A$ 的列数与 $B$ 的行数不同，则他们无法进行矩阵乘法。 是不是有点抽象？举个例子你就明白了。当 $A$ 的大小为 $2\times 2$ ， $B$ 的大小为 $2\times 3$ 时，则

$C = \begin{bmatrix} a_{11}b_{11}+a_{12}b_{21} & a_{11}b_{12}+a_{12}b_{22} & a_{11}b_{13}+a_{12}b_{23}\\ a_{21}b_{11}+a_{22}b_{21} & a_{21}b_{12}+a_{22}b_{22} & a_{21}b_{13}+a_{22}b_{23} \end{bmatrix} \text{.}$

同时，我们可以证明矩阵乘法满足结合律，但在一般情况下不满足交换律。

我们立刻写出了这样的矩阵乘法代码：

struct mat{
	int a[maxn][maxn];
};

mat mul(mat& a,mat& b){
	mat c;
	memset(c.a,0,sizeof(c.a));
	for(int k=1;k<=n;k++)
		for(int i=1;i<=n;i++)
			for(int j=1;j<=n;j++)
				c.a[i][j]+=a.a[i][k]*b.a[k][j];
	return c;
}

这个代码的时间复杂度是 $O(n^3)$ 的，可以接受。事实上，《算法导论》给出了一种时间复杂度 $O(n^{\log_2 7})$ 的矩阵乘法算法，但似乎并没有在 OI 界得到广泛运用，因此这里暂不介绍。

现在，我们回到正题：矩阵快速幂。由于矩阵乘法满足结合律，我们把它当作一般的快速幂打代码即可，只不过数字变成了矩阵，乘法变成了矩阵乘法。

实现

先放出 AC 模板题的代码：

#include <cstdio>
#include <algorithm>
#include <memory.h>

using namespace std;

inline long long read(){
	long long x=0,w=1;
	char ch=getchar();
	while(ch<'0'||ch>'9'){
		if(ch=='-') w=-1;
		ch=getchar();
	}
	while(ch>='0'&&ch<='9'){
		x=(x<<3)+(x<<1)+(ch^48);
		ch=getchar();
	}
	return x*w;
}

void write(int x){
	if(x<0) putchar('-'),x=-x;
	if(x>9) write(x/10);
	putchar((x%10)^48);
}

int n;

struct mat{
	int a[105][105];
};

mat mul(mat& a,mat& b){
	mat c;
	memset(c.a,0,sizeof(c.a));
	for(int k=1;k<=n;k++)
		for(int i=1;i<=n;i++)
			for(int j=1;j<=n;j++)
				c.a[i][j]=(c.a[i][j]+(0ll+a.a[i][k])*b.a[k][j]%1000000007)%1000000007;
	return c;
}

int main(){
	n=read();
	long long k=read();
	mat a;
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
			a.a[i][j]=read();
	
	mat ans;
	memset(ans.a,0,sizeof(ans.a));
	for(int i=1;i<=n;i++)ans.a[i][i]=1;
	
	do{
		if(k&1) ans=mul(ans,a);
		a=mul(a,a);k>>=1;
	}while(k);
	
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++){
			write(ans.a[i][j]);
			putchar(32);
		}
		putchar(10);
	}
	return 0;
}

整个代码的精髓是这几行：

mat ans;
memset(ans.a,0,sizeof(ans.a));
for(int i=1;i<=n;i++)ans.a[i][i]=1;

do{
    if(k&1) ans=mul(ans,a);
    a=mul(a,a);k>>=1;
}while(k);

最令人疑惑的是第三行。按照第三行的定义，我们构造出的矩阵是：

$Ans = \begin{bmatrix} 1 & 0 & \cdots & 0\\ 0 & 1 & \cdots & \vdots\\ \vdots & \vdots & \ddots & \vdots\\ 0 & \cdots & \cdots & 1\\ \end{bmatrix} \text{.}$

这个矩阵我们称其为单位矩阵。任何矩阵与单位矩阵相乘，得到的依然是原矩阵。 因此，单位矩阵就类似数的乘法中的 $1$ ，起到基的作用。

应用

例题：斐波那契数列（洛谷）或斐波那契数列（ACGO 搬运）。

题面非常简单，就是让你求斐波那契数列的第 $n$ 项。但是看到数据范围， $1\leq n<2^{63}$ ，显然事情并不是那么简单。现在，我们先从斐波那契数列的本质想起。它的递推式是这样的：

$f_i=f_{i-1}+f_{i-2}$

发现 $f_i$ 的值只与前两项有关。我们只需要保存两项的值就够了。我们把接下来需要保存的两项 $f_i,f_{i-1}$ 通过现在已知的两项 $f_{i-1},f_{i-2}$ 表达出来，形成一个线性方程组：

$\begin{cases} 1\times f_{i-1}+1\times f_{i-2} &= f_i\\ 1\times f_{i-1}+0\times f_{i-2} &= f_{i-1} \end{cases}$

根据矩阵乘法，我们可以把这个方程组表达成两个矩阵相乘：

$\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}f_{i-1}\\f_{i-2}\end{bmatrix}$ 也可以用这个式子求出来。于是，我们可以把它展开：

$\begin{bmatrix} 1 & 1\\ 1 & 0 \end{bmatrix}\left(\begin{bmatrix} 1 & 1\\ 1 & 0 \end{bmatrix}\begin{bmatrix} f_{i-2}\\ f_{i-3} \end{bmatrix}\right)=\begin{bmatrix} f_i\\ f_{i-1} \end{bmatrix}$

继续展开直到已知项 $\begin{bmatrix}f_2\\f_1\end{bmatrix}$ ：

$\begin{bmatrix} 1 & 1\\ 1 & 0 \end{bmatrix}\left(\begin{bmatrix} 1 & 1\\ 1 & 0 \end{bmatrix}\left(\cdots\begin{bmatrix} 1 & 1\\ 1 & 0 \end{bmatrix}\begin{bmatrix} f_2\\ f_1 \end{bmatrix}\cdots\right)\right)=\begin{bmatrix} f_i\\ f_{i-1} \end{bmatrix}$

又根据矩阵乘法的结合律：

$\begin{bmatrix} 1 & 1\\ 1 & 0 \end{bmatrix}^{i-2}\begin{bmatrix} f_2\\ f_1 \end{bmatrix}=\begin{bmatrix} f_i\\ f_{i-1} \end{bmatrix}$

（为什么是 $i-2$ 次？因为乘 $i$ 次的时候结果矩阵的首项就是 $f_{i+2}$ 了，所以乘 $i-2$ 次的时候结果矩阵的首项就刚好是 $f_i$ 了，如果不理解手推一下就明白了。）

现在，你一定发现了： $\begin{bmatrix}1 & 1\\1 & 0\end{bmatrix}^{i-2}$ 这一项，直接使用矩阵快速幂即可。最后再乘上 $\begin{bmatrix}f_2\\f_1\end{bmatrix}$ 即 $\begin{bmatrix}1\\1\end{bmatrix}$ ，得到的结果矩阵首项即为所求的 $f_i$ 。由于矩阵的大小为常数，所以可以把单次矩阵乘法视为 $O(1)$ ，矩阵快速幂的时间复杂度即为 $O(\log n)$ ，轻松跑过 long long 范围。注意特判 $n\leq 2$ 即可（本搬题人好心地加上了这两个极限情况）。AC 代码如下：

#include <cstdio>
#include <algorithm>
#include <memory.h>

using namespace std;

inline long long read(){
	long long x=0,w=1;
	char ch=getchar();
	while(ch<'0'||ch>'9'){
		if(ch=='-') w=-1;
		ch=getchar();
	}
	while(ch>='0'&&ch<='9'){
		x=(x<<3)+(x<<1)+(ch^48);
		ch=getchar();
	}
	return x*w;
}

void write(int x){
	if(x<0) putchar('-'),x=-x;
	if(x>9) write(x/10);
	putchar((x%10)^48);
}

struct mat{
	int a[3][3];
};

mat mul(mat& a,mat& b){
	mat c;
	memset(c.a,0,sizeof(c.a));
	for(int k=1;k<=2;k++)
		for(int i=1;i<=2;i++)
			for(int j=1;j<=2;j++)
				c.a[i][j]=(c.a[i][j]+(0ll+a.a[i][k])*b.a[k][j]%1000000007)%1000000007;
	return c;
}

int main(){
	long long n=read();
	if(n<=2){
		putchar(49);putchar(10);return 0;
	}
	n-=2;
	mat a;
	
	#define a a.a
	a[1][1]=1;a[1][2]=1;
	a[2][1]=1;a[2][2]=0;
	#undef a
	
	mat ans;
	memset(ans.a,0,sizeof(ans.a));
	ans.a[1][1]=1;ans.a[2][2]=1;
	
	do{
		if(n&1) ans=mul(ans,a);
		a=mul(a,a);n>>=1;
	}while(n);
	
	mat b;
	b.a[1][1]=b.a[2][1]=1;
	ans=mul(ans,b);
	
	write(ans.a[1][1]);putchar(10);
	return 0;
}

总结

总结一下，用矩阵快速幂优化此类递推（动态规划）的步骤：

确定哪些状态对结果状态或过程状态有贡献。只保留有贡献的状态。
确定上一次得到的状态转移到当前状态使用的递推式（方程组）。
提取方程组的各项系数为转移矩阵。
计算转移矩阵的次数。
写下矩阵快速幂和矩阵乘法模板。

接下来大家可以尝试一下矩阵加速（数列）（洛谷）、矩阵加速（数列）（ACGO 搬运）和广义斐波那契数列（洛谷）。如果你会数论，还可以尝试斐波那契公约数~~（我不会）~~。

祝大家暑假快乐！
~~看在我搬题，打 $\LaTeX{}$ ，写 Markdown 这么努力，可否留个赞支持一下 qwq。~~

Upd on 2024/6/30：ACGO主题库已有广义斐波那契数列与斐波那契公约数。

有帮助，赞一个

去预览

0/2000

全部评论 5

Macw07
算法导论提到的方法有点类似于分治的原理。

2024-06-22 来自浙江
1
- 复仇者_必胜啊
  回复Macw07
  矩阵快速幂……
  
  2024-07-27 来自广东
  0
Starsfocxy
顶

2024-08-23 来自江苏
0
暑假神（开学祭
顶

2024-06-24 来自上海
0
暑假神（开学祭
写完了：）

2024-06-23 来自上海
0
流萤是我的O.o
6
2024-06-23 来自河南
0

前置知识：

概念

实现

应用

总结

全部评论 5

热门讨论