矩阵乘法学习笔记

矩阵是一个数学阵列。

概念

一个 $m$ 行 $n$ 列的矩阵 $A$ 可以表示为：
$$
A=
egin{bmatrix}
a_{11} & a_{12} & cdots & a_{1n} \
a_{21} & a_{22} & cdots & a_{2n} \
cdots & cdots & & cdots \
a_{m1} & a_{m2} & cdots & a_{mn} \
end{bmatrix}
$$
这是一个 $m imes n$ 矩阵。

在这里，我们定义数组第一维为行号，第二维为列号，a[i][j] 表示矩阵第 $i$ 行第 $j$ 列的数。

矩阵乘法

当且仅当，第一个矩阵列数与第二个矩阵行数相等时，才可以将两者相乘。

设 $A$ 是 $m imes p$ 矩阵，$B$ 是 $p imes n$ 矩阵，则 $C = A imes B$ 是 $m imes n$ 矩阵。

其中 $C_{ij}$ 为 $A$ 第 $i$ 行第一个乘以 $B$ 第 $j$ 列第一个，加上 $A$ 第 $i$ 行第二个乘以 $B$ 第 $j$ 列第二个 …… 一直加到 $A$ 第 $i$ 行第 $p$ 个乘以 $B$ 第 $j$ 列第 $p$ 个。即：
$$
egin{aligned}
forall i in [1, m], j in [1, n] \
C_{ij} = sum_{k=1}^{p} A_{ik} imes B_{kj}
end{aligned}
$$
直接模拟即可。

性质

矩阵乘法满足结合律和分配律：
$$
egin{aligned}
(A imes B) imes C = A imes (B imes C)\
A imes C + B imes C = (A + B) imes C
end{aligned}
$$
但不满足交换律：
$$
A imes B ot= B imes A
$$

矩阵快速幂

根据矩阵乘法的性质，利用乘法结合律，基于快速幂可以快速求出 $A^n$。

模板

Luogu 3390

给定 $n imes n$ 的矩阵 $A$，求 $A^k$。

#include <cstdio>
#include <cstring>

typedef long long ll;

const int maxn = 105;
const int P = 1e9 + 7;

ll a[maxn][maxn];
ll t[maxn][maxn]; // t 为临时变量
ll ans[maxn][maxn];

void addToAns(int n) {
	memcpy(t, ans, sizeof(ans));
	memset(ans, 0, sizeof(ans));
	for (int i = 1; i <= n; i++) {
		for (int j = 1; j <= n; j++) {
			for (int k = 1; k <= n; k++) {
                // ans = ans * a
				ans[i][j] = (ans[i][j] + (t[i][k] * a[k][j]) % P) % P;
			}
		}
	}
}

void mulSelf(int n) {
	memcpy(t, a, sizeof(a));
	memset(a, 0, sizeof(a));
	for (int i = 1; i <= n; i++) {
		for (int j = 1; j <= n; j++) {
			for (int k = 1; k <= n; k++) {
                // a = a * a
				a[i][j] = (a[i][j] + (t[i][k] * t[k][j]) % P) % P;
			}
		}
	}
}

int main() {
	int n;
	ll k;
	scanf("%d %lld", &n, &k);
	for (int i = 1; i <= n; i++) {
		for (int j = 1; j <= n; j++) {
			scanf("%lld", &a[i][j]);
			ans[i][j] = a[i][j];
		}
	}
	k = k - 1;

	while (k) {
		if (k & 1) addToAns(n);
		mulSelf(n);
		k >>= 1;
	}

	for (int i = 1; i <= n; i++) {
		for (int j = 1; j <= n; j++) {
			printf("%lld ", ans[i][j]);
		}
		printf("
");
	}

	return 0;
}

例题

POJ 3070