[数据结构与算法-10]矩阵加速

矩阵加速

思路

利用矩阵运算，将递推数列转化为矩阵之幂，利用矩阵快速幂加速数列的计算

---

P1939 【模板】矩阵加速（数列）

题目描述

已知一个数列 a，它满足：

\[a_x = \left\{ \begin{matrix} 1 & x\in\{1, 2, 3\}\\ a_{x-1} + a_{x-3} & x \geq 4 \end{matrix} \right. \]

求 a 数列的第 n 项对 10^9+7 取余的值。

输入格式

第一行一个整数 T，表示询问个数。

以下 T 行，每行一个正整数 n。

输出格式

每行输出一个非负整数表示答案。

输入输出样例

输入 #1复制

3
6
8
10

输出 #1复制

4
9
19

说明/提示

• 对于 30% 的数据 n \(\leq\) 100 ；
• 对于 60% 的数据 n \(\leq2 \times 10^7\)；
• 对于 100% 的数据 \(1 \leq T \leq 100,1 \leq n \leq 2 \times 10^9\)。

解答

构造向量

\[\left(\begin{matrix} f[i]\\ f[i-1]\\ f[i-2] \end{matrix}\right) \]

寻找加速矩阵M

利用递推公式

\[\left(\begin{matrix} f[i] \\ f[i-1] \\ f[i-2] \\ \end{matrix}\right)= \left(\begin{matrix} f[i-1] \times 1 + f[i-2] \times 0 + f[i-3] \times 1\\ f[i-1] \times 1 + f[i-2] \times 0 + f[i-3] \times 0\\ f[i-1] \times 0 + f[i-2] \times 1 + f[i-3] \times 0\\ \end{matrix}\right)= \left(\begin{matrix} 1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{matrix}\right) \times \left(\begin{matrix} f[i-1] \\ f[i-2] \\ f[i-3] \\ \end{matrix}\right) \]

求出通项

\[\left(\begin{matrix} f[i] \\ f[i-1] \\ f[i-2] \\ \end{matrix}\right)= \left(\begin{matrix} 1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{matrix}\right)^{i - 1} \times \left(\begin{matrix} f[1] \\ f[2] \\ f[3] \\ \end{matrix}\right) \]

利用矩阵快速幂求解

#include <cstdio>
#include <cstring>
#define MOD 1000000007
int T, N;
struct Matrix {
	long long e[4][4];
}M,E;
// 重载运算符
Matrix operator*(const Matrix& a, const Matrix& b) {
	Matrix ans;
	memset(ans.e, 0, sizeof(ans.e));
	for (int i = 1; i <= 3; i++)
		for (int j = 1; j <= 3; j++)
			for (int k = 1; k <= 3; k++) {
				ans.e[i][j] += a.e[i][k] * b.e[k][j] % MOD;
				ans.e[i][j] %= MOD;
			}
	return ans;
}
int main() {
	scanf("%d", &T);
	for (int i = 0; i < T; i++) {
		scanf("%d", &N);
        // 初始化
		memset(E.e, 0, sizeof(E.e));
		memset(M.e, 0, sizeof(M.e));
		for (int i = 1; i <= 3; i++)
			E.e[i][i] = 1;
		M.e[1][1] = M.e[1][3] = M.e[2][1] = M.e[3][2] = 1;
        // 快速幂
		while (N) {
			if (N & 1)E = E * M;
			M = M * M;
			N >>= 1;
		}
		printf("%d\n", E.e[2][1]);
	}
	return 0;
}