矩阵快速幂

前置芝士：矩阵乘法，快速幂。

运算律

矩阵乘法不满足交换律，满足结合律。
所以用来加速递推的时候，我们就运用的是结合律。

以下公式 \(A\) 为 \(x \times y\) 矩阵，\(M\) 为构造出来的大小为 \(y \times z\) 矩阵，即：

\[\begin{split} A \times M \times M \times M \times \dots \times M &= A \times (M \times M \times M \times \dots \times M) \\ &= A \times M^K \end{split} \]

P3390 【模板】矩阵快速幂

这题是模版。

做法

矩阵求幂的话正常方法时间复杂度是 \(\mathcal O(x · y · z · k)\)

而快速幂加速乘法求幂的时间复杂度是 $\mathcal O(\log k) $。

那如果二者相结合，就可以将时间复杂度优化到 \(\mathcal O(x · y · z · \log k)\)

代码

戳我喵~

namespace MATRIX {
	const int mod = 1e9 + 7;
	template<typename Matrix>
	struct matrix {
		int size_n = 0, size_m = 0;
		int a[5][5];
		
		void resize(int n, int m) {
			size_n = n, size_m = m;
			for (int i = 0; i < size_n; i++) for (int j = 0; j < size_m; j++) a[i][j] = 0;
		}
		
		void init(int n) {
			size_n = n, size_m = n;
			for (int i = 0; i < size_n; i++) for (int j = 0; j < size_m; j++) a[i][j] = i == j ? 1 : 0;
		}
		
		auto& operator[](const int& i) {return a[i];}
		
		void input() {
			for (int i = 0; i < size_n; i++) for (int j = 0; j < size_m; j++) a[i][j] = re;
		}
		
		void output() {
			for (int i = 0; i < size_n; i++) {
				for (int j = 0; j < size_m; j++) wr((a[i][j] + mod) % mod), sp;
				endl;
			}
		}
		
		matrix friend operator*(matrix X, matrix Y) {
			assert(X.size_m == Y.size_n);
			matrix ans;
			ans.resize(X.size_n, Y.size_m);
			for (int i = 0; i < X.size_n; i++) for (int j = 0; j < Y.size_m; j++) for (int k = 0; k < X.size_m; k++) ans.a[i][j] = (ans.a[i][j] + 1ll * X.a[i][k] * Y.a[k][j] % mod) % mod;
			return ans;
		}
		
		friend matrix& operator*=(matrix &X, const matrix &Y) {
			X = X * Y;
			return X;
		}
		
		friend matrix operator%(matrix X, const Matrix &mod) {
			for (int i = 0; i < X.size_n; i++) for (int j = 0; j < X.size_m; j++) X[i][j] = (X[i][j] % mod + mod) % mod;
			return X;
		}
		
		friend matrix& operator%=(matrix& X, const Matrix &mod) {
			X = X % mod;
			return X;
		}
	};
	
	template<typename Matrix>
	
	matrix<Matrix> qpow(matrix<Matrix> a, int b) {
		matrix<Matrix> ans;
		ans.init(a.size_n);
		while (b) {
			if (b & 1) (ans *= a) %= mod;
			(a *= a) %= mod, b >>= 1;
		}
		return ans;
	}
}

---

那么再来一道模版吧！

P1962 斐波那契数列

不要高兴的太早， \(2 \le n \le 2^{63}\)

做法

我们推了一下可知一个矩阵 \(M\) 为：

\[M = \begin{bmatrix} 0 & 1 \\ 1 & 1 \end{bmatrix} \]

我们用一个矩阵 \(A\) 表示状态：

\[A = [f_n, f_{n-1}] \]

那么初始值就是

\[A = [1, 1] \]

由此可知，式子为：

\[A \times M^k \]

最后的答案就是 \(M_{0, 1}\)

为什么能省去一个乘 \(A\) 呢？就是因为 \(A = [1,1]\)

代码

戳我喵~

//上面就是模版
int main() {
    matrix<lint> M;
	M.resize(2, 2);
	M[0][1] = M[1][0] = M[1][1] = 1;
	int n = re;
	M = qpow(M, n);
	wr(M[0][1]), endl;
}

推出矩阵 \(M\) 的小妙招

构造步骤

确定需要维护的量：列出所有递推中需要用到的项

检查递推关系：确保这些量之间可以通过线性组合相互表示

设计状态向量：通常包含当前项、前几项、以及必要的辅助项

逐行构造矩阵：状态向量的每个分量对应矩阵的一行

验证：确保乘法后能得到正确的结果

常见技巧

常数项 \(\to\) 在状态向量中加入常数1

指数函数 \(c^n\) \(\to\) 加入 \(c^n\)，利用 \(c^{n} = c \cdot c^{n-1}\)

和式 \(S_n = \sum f(i)\) \(\to\) 加入 \(S_n\)，利用 \(S_n = S_{n-1} + f(n)\)

乘积项 \(a_n a_{n-1}\) \(\to\) 可能需要同时维护平方项和乘积项