矩阵乘法与矩阵快速幂

合集 - 学习笔记(52)

1.随机化算法2024-02-272.组合数学2024-02-273.概率和期望2024-02-274.树链剖分2024-02-27

5.矩阵乘法与矩阵快速幂2024-02-27

6.最短路2024-02-277.卡特兰数、Prüfer 序列、BSGS2024-02-278.分块和莫队2024-02-279.AC 自动机2024-02-2710.平衡树2024-02-2711.基础字符串算法2024-02-2712.基础数论2024-02-2713.Miller-Rabin 与 Pollard-Rho2024-08-2314.广义后缀自动机2024-06-2515.后缀自动机 SAM2024-06-2216.回文自动机 PAM2024-06-1817.Manacher2024-06-1618.后缀数组2024-06-0819.01 分数规划2024-06-0220.网络流2024-05-0521.整体二分2024-05-0522.cdq 分治2024-05-0423.点分治2024-04-2124.虚树2024-04-1425.扫描线2024-04-1126.笛卡尔树2024-04-0627.基环树2024-04-0228.树哈希2024-03-0929.莫比乌斯反演2024-11-2830.二分图2024-11-0731.朱刘算法2024-10-2932.最大团2024-10-2933.杂项2024-08-2434.杜教筛2024-12-0135.拉格朗日插值2024-12-1036.线段树综合2024-12-1737.可持久化数据结构2024-12-2438.K-D Tree2024-12-2539.Burnside 引理与 Polya 定理2024-12-2740.线性基2024-12-3041.替罪羊树01-0942.LCT01-1843.插头 dp01-2044.原根01-2045.多项式乘法01-2246.斯特林数01-2347.二项式反演与斯特林反演01-2348.Min-Max 容斥01-2449.辛普森积分法01-2550.Min_25 筛02-0351.凸包02-1652.2-SAT03-03

收起

1 矩阵乘法

1.1 定义

对于两个矩阵 $A,B$，其中 $A$ 大小为 $n\times m$ ，$B$ 大小为 $m\times p$，则这两个矩阵可以做乘法，得到的矩阵 $C$ 的大小为 $n\times p$。

例如：

$A=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\end{array}\right]$

$B=\left[\begin{array}{cc}{b}_{11}& b_{12}\\ b_{21}& b_{22}\\ b_{31}& b_{32}\end{array}\right]$

则有：

$$C=AB=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\end{array}\right]\left[\begin{array}{cc}{b}_{11}& b_{12}\\ b_{21}& b_{22}\\ b_{31}& b_{32}\end{array}\right]=\left[\begin{array}{cc}{a}_{11}{b}_{11}+a_{12}{b}_{21}+a_{13}{b}_{31}& a_{11}{b}_{12}+a_{12}{b}_{22}+a_{13}{b}_{32}\\ a_{21}{b}_{11}+a_{22}{b}_{21}+a_{23}{b}_{31}& a_{21}{b}_{12}+a_{22}{b}_{22}+a_{23}{b}_{32}\end{array}\right]$$

1.2 注意点

• 只有当矩阵 $A,B$ 满足 $A$ 的列数等于 $B$ 的行数时，才能进行矩阵乘法。
• 矩阵乘法不满足交换律（这点很重要）

1.3 代码

提供一个封装模板。

#include <bits/stdc++.h>
#define int long long

using namespace std;

typedef long long LL;
const int Maxn = 205;

int n, m, p;

struct matrix {
	int n, m, p[Maxn][Maxn];
	matrix operator * (const matrix &b) const {
		matrix a = *this, c;
		c.n = a.n, c.m = b.m;
		for(int i = 1; i <= a.n; i++) {
			for(int j = 1; j <= a.m; j++) {
				for(int k = 1; k <= b.m; k++) {
					c.p[i][k] += a.p[i][j] * b.p[j][k];
				}
			}
		}
		return c;
	}
	void print() {
		for(int i = 1; i <= n; i++) {
			for(int j = 1; j <= m; j++) {
				cout << p[i][j] << " ";
			}
			cout << '\n';
		}
	}
}A, B, C;

signed main() {
	ios::sync_with_stdio(0);
	cin >> n >> m;
	A.n = n, A.m = m;
	for(int i = 1; i <= n; i++) {
		for(int j = 1; j <= m; j++) {
			cin >> A.p[i][j];
		}
	}
	cin >> p;
	B.n = m, B.m = p;
	for(int i = 1; i <= m; i++) {
		for(int j = 1; j <= p; j++) {
			cin >> B.p[i][j];
		}
	}
	C = A * B;
	C.print();
	return 0;
}

2 矩阵快速幂

2.1 定义

仿照着数字的快速幂，我们也能快速求出矩阵的幂。

定义矩阵的幂为：$A^n=A\cdot A\cdot A\cdot A\cdots$

这很明显需要满足原矩阵为正方形。

2.2 用途

矩阵快速幂的经典应用就是加速递推。

直接举几个例子。

2.2.1 「例 1」Fibonacci 第 n 项

Problem

求斐波那契数列第 $n$ 项后四位（$n\le {10}^9$）

Solution

第一眼：这不是递推板题吗！

然而 $n\le {10}^9$，$O(n)$ 暴力绝对爆炸。

我们运用矩阵快速幂加速递推。

接下来我们来求一下递推式（注意掌握方法）：

$\mathrm{①}$ 我们先来写出 $F_n$ 的部分：

$$\left[\begin{array}{c}{F}_n\\ \end{array}\right]=\left[\begin{array}{cc}1& 1\\ & \end{array}\right]\left[\begin{array}{c}{F}_{n-1}\\ F_{n-2}\end{array}\right]$$

$\mathrm{②}$ 接下来我们关注右边的第二个矩阵，发现他们的下标相差 $1$，因此左边的矩阵应该是 $\left[\begin{array}{c}{F}_n\\ F_{n-1}\end{array}\right]$ 。

$\mathrm{③}$ 最后补全整个式子：

$$\left[\begin{array}{c}{F}_n\\ F_{n-1}\end{array}\right]=\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]\left[\begin{array}{c}{F}_{n-1}\\ F_{n-2}\end{array}\right]$$

因此我们对于矩阵 $\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]$ 不断做快速幂，然后乘上 $\left[\begin{array}{c}{F}_1\\ F_0\end{array}\right]$ 即可。

Code

#include <bits/stdc++.h>
#define int long long

using namespace std;

typedef long long LL;
const int Maxn = 205;

int n, mod;

struct matrix {
	int n, m, p[Maxn][Maxn];
	void init() {
		n = 0, m = 0, memset(p, 0, sizeof p);
	}
	matrix operator * (const matrix &b) const {
		matrix a = *this, c;
		c.init();
		c.n = a.n, c.m = b.m;
		for(int i = 1; i <= a.n; i++) {
			for(int j = 1; j <= a.m; j++) {
				for(int k = 1; k <= b.m; k++) {
					c.p[i][k] = (c.p[i][k] + a.p[i][j] * b.p[j][k] % mod) % mod;
				}
			}
		}
		return c;
	}
	matrix operator ^ (const int &o) const {
		matrix res, a = *this; int b = o;
		res.init();
		res.n = res.m = 2;
		res.p[1][1] = res.p[2][2] = 1;
		while(b) {
			if(b & 1) {
				res = res * a;
			}
			a = a * a;
			b >>= 1;
		}
		return res;
	}
	void print() {
		for(int i = 1; i <= n; i++) {
			for(int j = 1; j <= m; j++) {
				cout << p[i][j] << " ";
			}
			cout << '\n';
		}
	}
}A;

signed main() {
	ios::sync_with_stdio(0);
	cin >> n >> mod; 
	A.p[1][1] = A.p[1][2] = A.p[2][1] = 1;
	A.n = A.m = 2;
	A = A ^ (n - 1);
	cout << A.p[1][1];
	return 0;
}

2.2.2 「例 2」P207 迷路

Solution

我们发现如果边权为 $1$ ，那么直接将矩阵 $T$ 次方即可。

但是此时边权不为 $1$，所以不能这样做。

但是我们又会发现这个边权居然最高才 $9$，因此我们考虑一下分层图的思路。

将一个点拆成 $9$ 个点并顺次相连，这样能将边权转化为 $1$。举个例子：

对于这样一条边，我们将 $1$ 和 $2$ 两个点拆开，然后将他们按照原先边权连起来，如下图：

这样从 $1_1$ 走到 $2_1$，经过的边权仍然是 $6$，但我们就转化成了边权为 $1$ 的图。

注意最后统计答案时统计的是点 $n_1$ 的值而非 $n$ 的所有值的和。

接下来注意细节即可。

Code

#include <bits/stdc++.h>

/*
我们发现如果边权为 1，那就是一道求 t 次方的大板子题
然而这题有边权，所以不能这么干
我们考虑分层图的思想，既然边权才 9，那我们直接将一个点暴力拆成九个点，然后边权就能转化为 1
这下就是大板子题了，在求 t 次方即可 
*/

using namespace std;

typedef long long LL;
const int Maxn = 205;
const int Mod = 2009;

int n, t, m;

struct matrix {
	int n, m, p[Maxn][Maxn];
	void init() {
		n = 0, m = 0, memset(p, 0, sizeof p);
	}
	matrix operator * (const matrix &b) const {
		matrix a = *this, c;
		c.init();
		c.n = a.n, c.m = b.m;
		for(int i = 1; i <= a.n; i++) {
			for(int j = 1; j <= a.m; j++) {
				for(int k = 1; k <= b.m; k++) {
					c.p[i][k] = (c.p[i][k] + a.p[i][j] * b.p[j][k] % Mod) % Mod;
				}
			}
		}
		return c;
	}
	matrix operator ^ (const int &o) const {
		matrix res, a = *this; int b = o;
		res.init();
		res.n = res.m = a.n;
		for(int i = 1; i <= n; i++) {
			res.p[i][i] = 1;
		}
		while(b) {
			if(b & 1) {
				res = res * a;
			}
			a = a * a;
			b >>= 1;
		}
		return res;
	}
	void print() {
		for(int i = 1; i <= n; i++) {
			for(int j = 1; j <= m; j++) {
				cout << p[i][j] << " ";
			}
			cout << '\n';
		}
	}
}A;

signed main() {
	ios::sync_with_stdio(0);
	cin >> n >> t;
	m = n * 9;
	for(int i = 1; i <= n; i++) {
		for(int j = 1; j <= 8; j++) {
			A.p[9 * (i - 1) + j][9 * (i - 1) + j + 1] = 1;
		}
	}
	for(int i = 1; i <= n; i++) {
		for(int j = 1; j <= n; j++) {
			char ch;
			cin >> ch;
			if(ch > '0') {
				A.p[9 * (i - 1) + ch - '0'][9 * (j - 1) + 1] = 1;
			}
		}
	}
	A.n = A.m = m;
	A = A ^ t;
	cout << A.p[1][9 * (n - 1) + 1];
	return 0;
}

2.2.3 「例 3」佳佳的 Fibonacci

Problem

求 $(F_1+2F_2+3F_3+\cdots +n{F}_n)modm$ 的值。

Solution

我们还是运用矩阵快速幂。

但是我们发现这个形式似曾相识，还记得求区间修改区间查询的 BIT 时候的式子吗？他和这个形式是一样的。

我们还是再来推一遍：

先设 $S_i=F_1+F_2+\cdots +F_i$

然后有：

$$\begin{array}{rl}& F_1+2F_2+3F_3+\cdots +n{F}_n\\ & =n\times S_n-S_{n-1}-S_{n-2}-\cdots -S_1\\ & =n\times S_n-\sum _{i=1}^{n-1}{S}_i\end{array}$$

然后我们再令 $P_n=\sum _{i=1}^{n-1}{S}_i$，则此时原式就可以化为：

$$n\times S_n-P_n$$

接下来我们递推 $S_n,P_n$ 来得到原式。

还是分步讨论：

$\mathrm{①}$ 写出关于 $P_n$ 的部分。注意到 $P_n=P_{n-1}+S_{n-1}$：

$\left[\begin{array}{c}{P}_n\\ \\ \\ \end{array}\right]=\left[\begin{array}{cccc}1& 1& 0& 0\\ & & & \\ & & & \\ & & & \end{array}\right]\left[\begin{array}{c}{P}_{n-1}\\ S_{n-1}\\ \\ \end{array}\right]$

$\mathrm{②}$ 注意到右边第二个矩阵要满足 $S$ 与 $P$ 下标相等，所以左边还要满足。又注意到 $S_n=S_{n-1}+F_n$，所以再补全一部分：

$\left[\begin{array}{c}{P}_n\\ S_n\\ \\ \end{array}\right]=\left[\begin{array}{cccc}1& 1& & \\ 0& 1& 1& \\ & & & \\ & & & \end{array}\right]\left[\begin{array}{c}{P}_{n-1}\\ S_{n-1}\\ F_n\\ \end{array}\right]$

$\mathrm{③}$ 又注意到右边 $F$ 的下标比 $S$ 的下标多 $1$，因此左边应该是 $F_{n+1}$。又注意到 $F_{n+1}=F_n+F_{n-1}$，所以补全为：

$\left[\begin{array}{c}{P}_n\\ S_n\\ F_{n+1}\\ F_n\end{array}\right]=\left[\begin{array}{cccc}1& 1& 0& 0\\ 0& 1& 1& 0\\ 0& 0& 1& 1\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{c}{P}_{n-1}\\ S_{n-1}\\ F_n\\ F_{n-1}\end{array}\right]$

然后实现即可。

Code

十年 OI 一场空，不开 long long 见祖宗。

#include <bits/stdc++.h>
#define int long long

using namespace std;

typedef long long LL;
const int Maxn = 205;

int n, mod;

struct matrix {
	int n, m, p[Maxn][Maxn];
	void init() {
		n = m = 0;
		memset(p, 0, sizeof p);
	}
	void print() {
		for(int i = 1; i <= n; i++) {
			for(int j = 1; j <= m; j++) {
				cout << p[i][j] << " ";
			}
			cout << '\n';
		}
	}
	matrix operator * (const matrix &b) const {
		matrix a = *this, c;
		c.init();
		c.n = a.n, c.m = b.m;
		for(int i = 1; i <= a.n; i++) {
			for(int j = 1; j <= a.m; j++) {
				for(int k = 1; k <= b.m; k++) {
					c.p[i][k] = (c.p[i][k] + a.p[i][j] * b.p[j][k] % mod) % mod;
				}
			}
		}
		return c;
	}
	matrix operator ^ (const int &b) const {
		int p = b;
		matrix a = *this, res;
		res.init();
		res.n = res.m = a.n;
		for(int i = 1; i <= a.n; i++) {
			res.p[i][i] = 1;
		}
		while(p) {
			if(p & 1) {
				res = res * a;
			}
			a = a * a;
			p >>= 1;
		} 
		return res;
	}
}A, F;

signed main() {
	ios::sync_with_stdio(0);
	cin >> n >> mod;
	A.init(), F.init();
	A.n = A.m = 4;
	A.p[1][1] = A.p[1][2] = A.p[2][2] = A.p[2][3] = A.p[3][3] = A.p[3][4] = A.p[4][3] = 1;
	F.n = 4, F.m = 1;
	F.p[3][1] = 1;
	A = A ^ n;
	A = A * F;
	cout << (n * A.p[2][1] % mod - A.p[1][1] + mod)	% mod;
	return 0;
}