浅谈矩阵快速幂

矩阵快速幂

一、矩阵乘法

1. 定义

设 \(A=\left[a_{i,j}\right]_{m\times s}\)，\(B=\left[b_{i,j}\right]_{s\times n}\)，则

\[c_{i,j}=\sum^s_{k=1}a_{i,k}\times b_{k, j} \]

2. 运算律

1. 满足结合律，但不满足交换律；
2. 左矩阵列数=右矩阵行数，乘法才有意义；
3. E相当于数字 \(1\) 在整数乘法中的作用，设为 \(\left[e_{i,j}\right]_{s\times s}\)，则：

\[\begin{equation} E=\left[e_{i,j}\right]_{s\times s}= \begin{cases} 1& \text{ $ i=j $} \\ 0& \text{ $ i\not=j $} \end{cases} \end{equation} \]

4. 矩阵A和它同阶的单位矩阵E作乘积，结果仍为A，即 \(AE=EA=A\)。

3. 代码

#include <stdio.h>
#include <string.h>
#include <math.h>
const int N = (int)1e4 + 5;
const double pi = acos(-1.0);
const int SZ = 5;
struct Matrix {
	int m[SZ + 1][SZ + 1];
	inline void init() {
		memset(m, 0, sizeof m);
	}
	static Matrix E() {
		Matrix e; e.init();
		for (int i = 1; i <= SZ; ++i)
			for (int j = 1; j <= SZ; ++j)
				e.m[i][j] = (i == j);
		return e;
	}
	friend Matrix operator *(const Matrix &a, const Matrix &b) {
		Matrix Res; Res.init();
		for (int i = 1; i <= SZ; ++i)
			for (int j = 1; j <= SZ; ++j)
				for (int k = 1; k <= SZ; ++k)
					Res.m[i][j] += a.m[i][k] * b.m[k][j];
		return Res;
	}
};
int main(void) {
	Matrix e; e = Matrix:: E(); e = e * e;
	for (int i = 1; i <= 5; ++i) {
		for (int j = 1; j <= 5; ++j)
			printf("%d ", e.m[i][j]);
		putchar('\n');
	}
	return 0;
}

二、矩阵快速幂

1. 代码

int res[N][N];
inline void Pow(int a[][N], int n) {
	memset(res, 0, sizeof res);
	for (int i = 1; i <= n; ++i) res[i][i] = 1;
	while (n) {
		if (n & 1) multi(res, a, n);
		multi(a, a, n);
		n >>= 1;
	}
}

三、应用举例

例1、矩阵幂级数

分析：\(S(k)=A+S(k-1)\times A\)。
设 \(B\) 矩阵满足 \(\left[A\ S(k-1)\right]\times B=\left[A\ S(k)\right]\)，则有：

	\(A\)	\(S_k\)
\(A\)	\(E\)	\(E\)
\(S_{k-1}\)	\(0\)	\(A\)

例2、灯环的状态 HDU2276

分析：矩阵加法转为矩阵异或。

	\(1\)	\(2\)	\(3\)	\(\cdots\)	\(n\)
\(1\)	\(1\)	\(1\)	\(0\)	\(\cdots\)	\(0\)
\(2\)	\(0\)	\(1\)	\(1\)	\(\cdots\)	\(0\)
\(3\)	\(0\)	\(0\)	\(1\)	\(\cdots\)	\(0\)
\(\cdots\)	\(0\)	\(0\)	\(0\)	\(\cdots\)	\(0\)
\(n-1\)	\(0\)	\(0\)	\(0\)	\(\cdots\)	\(1\)
\(n\)	\(1\)	\(0\)	\(0\)	\(\cdots\)	\(1\)

例3、平面上的点

分析：按操作构造以下矩阵并作叉乘。

1. 平移 \((p, q)\)

	\(1\)	\(x\)	\(y\)
\(1\)	\(1\)	\(p\)	\(q\)
\(x\)	\(0\)	\(1\)	\(0\)
\(y\)	\(0\)	\(0\)	\(1\)

2. 缩放（\(L\)）

	\(1\)	\(x\)	\(y\)
\(1\)	\(1\)	\(0\)	\(0\)
\(x\)	\(0\)	\(L\)	\(0\)
\(y\)	\(0\)	\(0\)	\(L\)

3. 上下翻转

	\(1\)	\(x\)	\(y\)
\(1\)	\(1\)	\(0\)	\(0\)
\(x\)	\(0\)	\(1\)	\(0\)
\(y\)	\(0\)	\(0\)	\(-1\)

4. 左右翻转

	\(1\)	\(x\)	\(y\)
\(1\)	\(1\)	\(0\)	\(0\)
\(x\)	\(0\)	\(-1\)	\(0\)
\(y\)	\(0\)	\(0\)	\(1\)

5. 旋转 \(\alpha\)

	\(1\)	\(x\)	\(y\)
\(1\)	\(1\)	\(0\)	\(0\)
\(x\)	\(0\)	\(\cos\alpha\)	\(\sin\alpha\)
\(y\)	\(0\)	\(\sin\alpha\)	\(\cos\alpha\)