矩阵快速幂、矩阵加速递推

更加洛谷的观感感受
写数论文章写爽了。

矩阵快速幂

矩阵乘法

矩阵乘法要求

设 $A$ 为 $m \times p$ 的矩阵，$B$ 为 $p \times n$ 的矩阵。
矩阵乘法必须第一个矩阵的列数等于第二个矩阵的行数。
$A \times B$ 是一个 $m \times n$ 的矩阵。

矩阵乘法定义

设 $A \times B = C$，则 $C_{i,j}$ 为：

$$C_{i,j} = \sum_{k=1}^p A_{i,k} \times B_{k,j}$$

矩阵乘法示例

假如：

$$A = \begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \end{bmatrix}, \quad B = \begin{bmatrix} 1 & 2 \\ 3 & 4 \\ 5 & 6 \end{bmatrix}$$

请你计算一下 $A \times B$ 的值。

我们一起计算一下：

$$\begin{aligned} C_{1,1} &= 1 \times 1 + 2 \times 3 + 3 \times 5 = 22 \\ C_{1,2} &= 1 \times 2 + 2 \times 4 + 3 \times 6 = 28 \\ C_{2,1} &= 4 \times 1 + 5 \times 3 + 6 \times 5 = 49 \\ C_{2,2} &= 4 \times 2 + 5 \times 4 + 6 \times 6 = 64 \end{aligned}$$

即：

$$C = \begin{bmatrix} 22 & 28 \\ 49 & 64 \end{bmatrix}$$

矩阵的幂运算

计算 $C = \begin{bmatrix}1 & 2 \\ 2 & 1\end{bmatrix}^3$：

1. 先算二次幂：

   $$\begin{bmatrix}1 & 2 \\ 2 & 1\end{bmatrix}^2 = \begin{bmatrix}1 & 2 \\ 2 & 1\end{bmatrix} \times \begin{bmatrix}1 & 2 \\ 2 & 1\end{bmatrix} = \begin{bmatrix}5 & 4 \\ 4 & 5\end{bmatrix}$$

2. 再算三次幂：

   $$\begin{bmatrix}5 & 4 \\ 4 & 5\end{bmatrix} \times \begin{bmatrix}1 & 2 \\ 2 & 1\end{bmatrix} = \begin{bmatrix}13 & 14 \\ 14 & 13\end{bmatrix}$$

单位矩阵

单位矩阵 $I = \begin{bmatrix}1 & 0 \\ 0 & 1\end{bmatrix}$，与任意矩阵相乘不改变其值。

一个 $n \times n$ 的单位矩阵为主对角线为 $1$，其他为 $0$ 的矩阵。

矩阵快速幂的实现

比较简单，建立一个矩阵结构体封装即可。

class matrix{
public:
	matrix()
		:n(0)
	{
		memset(c, 0, sizeof(c));
	}
	int* operator[](int t)
	{
		return c[t];
	}
	void init()//建立单位矩阵
	{
		for(int i = 1; i <= n; i++)
			c[i][i] = 1;
	}
	int size()
	{
		return n;
	}
private:
	int c[110][110];
	int n;
};
matrix operator*(matrix x, matrix y)//矩阵乘法
{
	matrix ans;
	int n = x.size();
	for(int i = 1; i <= n; i++)
		for(int j = 1; j <= n; j++)
			for(int k = 1; k <= n; k++)
				ans[i][j] += x[i][k] * y[k][j];
}
matrix quick_pow(matrix x, int y)//矩阵快速幂
{
	matrix ans;
	ans.init();//注意这里需要是单位矩阵
	while(y)
	{
		if(y & 1) ans = ans * x;
		y >>= 1;
		x = x * x;
	}
	return ans;
}

时间复杂度为：$O(n^3\times\log n )$，其中 $y$ 为指数， $n$ 为矩阵长度。

矩阵加速递推

求斐波那契数列（P1962 斐波那契数列）：

$$F_n=\begin{cases}1 & (n \leq2)\\ F_{n-1}+F_{n-2} & (n > 2)\end{cases}$$

$$1\leq n \leq2^{63}$$

思路

构造矩阵

我们把斐波那契数列的相邻两项表示为一个矩阵 $\begin{bmatrix}F_n & F_{n - 1}\end{bmatrix}$，我们希望通过 $\begin{bmatrix}F_{n - 1} & F_{n - 2}\end{bmatrix}$ 推出 $\begin{bmatrix}F_n & F_{n - 1}\end{bmatrix}$。

构造转移矩阵

我们构造一个矩阵 $A$，使得：

$$\begin{bmatrix}F_{n - 1} & F_{n - 2}\end{bmatrix}\times A=\begin{bmatrix}F_n & F_{n - 1}\end{bmatrix}$$

根据递推关系：

$$\begin{cases}F_n=F_{n-1}\times 1+F_{n-2}\times 1\\ F_{n-1}=F_{n-1}\times 1+F_{n-2}\times 0\end{cases}$$

观察可得：

$$A=\begin{bmatrix}1 & 1\\1 & 0\end{bmatrix}$$

递推

我们发现：

$$\begin{bmatrix}F_{2} & F_{1}\end{bmatrix}\times A^{n-2}=\begin{bmatrix}F_n & F_{n - 1}\end{bmatrix}$$

时间复杂度

矩阵快速幂为 $O(n^3 \log y)$，而矩阵的长度为 $2$，$y$ 为题目中的 $n - 2$，所以时间复杂度为 $O(8\log n)$，就是 $O(\log n)$。

代码

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int mod = 1e9 + 7;
int n = 2;
struct matrix{
	int c[3][3];
	matrix()
	{
		memset(c, 0, sizeof(c));
	}
	void init()
	{
		for(int i = 1; i <= n; i++)
			c[i][i] = 1;
	}
};
matrix operator*(matrix x, matrix y)
{
	matrix ans;
	for(int i = 1; i <= n; i++)
		for(int j = 1; j <= n; j++)
			for(int k = 1; k <= n; k++)
				ans.c[i][j] = (ans.c[i][j] + x.c[i][k] * y.c[k][j]) % mod;
	return ans;
}
matrix quick_pow(matrix x, int y)
{
	matrix ans;
	ans.init();
	while(y)
	{
		if(y & 1) ans = ans * x;
		y >>= 1;
		x = x * x;
	}
	return ans;
}
signed main()
{
	int m;
    cin >> m;
    if(m <= 2) cout << 1, exit(0);
    matrix a;
    a.c[1][1] = a.c[1][2] = 1;
    matrix b;
    b.c[1][1] = b.c[1][2] = b.c[2][1] = 1;
    b = quick_pow(b, m - 2);
    matrix ans = a * b;
    cout << ans.c[1][1];
	return 0;
}

作业

已知一个函数为：

$$f(n)=\begin{cases}1 &(n \leq 3)\\ f(n - 1) \times f(n - 2) \times f(n-3) & (n > 3)\end{cases}$$

给出 $n(1\leq n\leq10^{18})$ 求出 $f(n)\mod 114514$。