矩阵快速幂求解递推式问题进阶(超详细)

最新推荐文章于 2024-02-25 23:46:06 发布

Cookie_coolkid

最新推荐文章于 2024-02-25 23:46:06 发布

阅读量302

点赞数 6

分类专栏：算法文章标签：矩阵算法线性代数

本文链接：https://blog.csdn.net/weixin_69763109/article/details/129839831

版权

算法专栏收录该内容

15 篇文章 1 订阅

订阅专栏

简介：

在上一篇的blog中，我们在求解斐波那契数列问题上发现了新方法，即使是很大的n值也可以十分迅速求解其在模998244353空间中的数值。因为矩阵快速幂的算法将其运算时间复杂度降低到了 $\mathbf{O}(lg n)$ .（尽管在上一篇的blog中写出的代码并不简洁美观。在接触这个算法的思想之后，我尝试自己实现初始化矩阵以及实现矩阵快速幂的函数，十分的冗长，并且不能够保证其完全正确性。）

在本篇blog中会继续介绍矩阵快速幂算法的运用。分别是来自两道题，一道是是求解更加复杂的递推式（不仅包含线性递推项，还包括常数的质数次幂以及多项式表达式）的数值问题，还有一道题目是从所谓一维范围内的递推式拓展到了二维的递推式。来自于宽搜题马走日的改版（在前面的一篇博客中有提到），变成了一道完全不一样的题目(在本篇博客中仅展示代码，思路与上一题类似，只需要考虑如果进行二维到一维的转化即可)。

需要注意的是：取模也是一个很大问题，我在这次问题中少取一次多余的取模操作，时间竟然可以优化600ms（在某个比较复杂的输入情况下）可见我们也需要一个简洁有效的取模方式。

题目描述：

在这个问题下，我们的递推式变得更加复杂了。可以说矩阵快速幂函数是一个模板，但是初始化过度矩阵的过程才是最为关键的一步。我们可以分开来处理：

对于线性递推项部分：类比求解斐波那契的思路，在这里只是递推式的数量可能大于斐波那契数列，但本质的思想并没有改变。

例如递推式： $\large f_n = a_1f_{n-1} + a_2f_{n-2} +a_3f_{n-3}$

我们给出的初始矩阵是 $\large \begin{pmatrix} f_2 & f_1 & f_0 \end{pmatrix}$

对于 $\large \begin{pmatrix} f_{n-1}\\f_{n-2} \\ f_{n-3} \end{pmatrix}$ 我们希望在一次矩阵乘法后可以得到 $\large \begin{pmatrix} f_n\\f_{n-1} \\ f_{n-2} \end{pmatrix}$ 因此我们主要需要构造的是求解 $f_n$ 而剩下的只要是令他等于其前一项即可。（就像在设计贪吃蛇游戏一样，你只需要考虑头部下一格往哪里走，接下来的身体部分就需要不断向前赋值即可）因此我们可以构造这样的矩阵（对初始矩阵做右乘）：

$\large \begin{pmatrix} a_1& a_2 &a_3 \\ 1 &0 &0 \\ 0 &1 &0 \end{pmatrix}$

对于其中常数的指数次幂 $\large c^{n-1}$ 我们希望在一次矩阵乘法过后可以得到 $\large c^n$

例如在初始矩阵中 $\large \begin{pmatrix} 0\\ c \\ 0 \end{pmatrix}$ 我们要得到 $\large \begin{pmatrix} 0\\ c^2 \\ 0 \end{pmatrix}$ 只需要乘上这个矩阵即可： $\large \begin{pmatrix} 0 & 0& 0\\ 0& c& 0\\ 0&0 & 0 \end{pmatrix}$ (该矩阵在左，右乘上初始矩阵)

对于多项式表达式，例如 $\large \begin{pmatrix} b_1\\b_2n \\ b_3n^2 \end{pmatrix}$ 我们考虑初始矩阵为 $\large \begin{pmatrix} b_1\\ n \\ n^2 \end{pmatrix}$ 在我们构造矩阵的时候会遇到一些问题，在于这边的n的二次方项上：我们会发现若想从 $\large n^2$ 得到 $\large (n+1)^2$ 中间相差了 $\large 2n + 1$ 由于还需要考虑常数 $\large b_3$ 的存在，我们必须要补上两项 $\large b_3n$ 以及 $\large b_3$ 同理对于 n 的一次方项我们还需要补上 $\large b_2$ .因此出于这样考虑我们必须改变初始矩阵变成这样：

$\large \begin{pmatrix} b_1\\b_2 \\ b_2n \\ b_3 \\ b_3n \\ b_3n^2 \end{pmatrix}$

我们同样可以根据其系数特征来构造这样的矩阵：

$\large \begin{pmatrix} 1 &0 & 0&0 & 0 &0 \\ 0& 1 & 0 & 0& 0&0 \\ 0& 1 & 1& 0 & 0 &0 \\ 0& 0& 0 &1 & 0 &0 \\ 0& 0& 0& 1 &1 &0 \\ 0& 0& 0 & 1& 2 & 1 \end{pmatrix}$

到这里，我们就已经处理完了这边的所有情况啦~接下来就是把根据输入在这些合并成一个过渡矩阵既可，同时，这边给出矩阵乘法和矩阵快速幂的简洁写法（在这里我们推荐使用类来构造这样一个矩阵所拥有的方法，变得更加简洁）：

#define MOD 998244353
typedef unsigned long long ull;

class Mat {
public:
	ull matrix[30][30];
};

Mat init, out;

Mat Mul(Mat x, Mat y)
{
	Mat c;
	for (int i = 1; i <= all; ++i)
		for (int j = 1; j <= all; ++j)
			c.matrix[i][j] = 0;

	for (int i = 1; i <= all; ++i)
		for (int j = 1; j <= all; ++j)
			for (int k = 1; k <= all; ++k)
				c.matrix[i][j] = c.matrix[i][j] % MOD + x.matrix[i][k] * y.matrix[k][j] % MOD;
	return c;
}

Mat pow(Mat x, ull y)//矩阵快速幂
{
	Mat ans = init;
	while (y) {
		if (y & 1) ans = Mul(ans, x);
		x = Mul(x, x);
		y >>= 1;
	}
	return ans;
}

构造过程大家可以自己去尝试尝试~

题目描述：

这边仅给出参考代码（关键在于如何将二维转化为一维的递推式进行求解）：

#include<iostream>
using namespace std;

#define MOD 998244353
typedef long long ll;

const int M = 250;

class Mat
{
public:
	ll(* grid)[M] = new ll[M][M];
};

int m, n;
int x, y;
int k;
int dx[9] = { 0, -1, -1, 1, 1, -2, -2, 2, 2 };
int dy[9] = { 0, -2, 2, -2, 2, -1, 1,-1, 1 };
ll arr[M];
ll ans[M];
Mat init;

Mat mul(Mat x, Mat y)
{
	Mat c;
	for (int i = 1; i <= n * m; i++)
		for (int j = 1; j <= n * m; j++)
			c.grid[i][j] = 0;
	for (int i = 1; i <= n * m; i++)
	{
		for (int k = 1; k <= n * m; k++)
		{
			int j;
			for (j = 1; j <= n * m; j++)
				c.grid[i][j] = (c.grid[i][j] + x.grid[i][k] * y.grid[k][j]) % MOD;
		}
	}
	return c;
}

Mat qpow(Mat x, int n)
{
	Mat tmp = init;
	while (n)
	{
		if (n & 1) tmp = mul(tmp, x);
		x = mul(x, x);
		n >>= 1;
	}
	return tmp;
}

bool check(int x, int y, int dx,int dy)
{
	if (x + dx < 1 || x + dx > m || y + dy < 1 || y + dy > n)
		return false;
	return true;
}

void initialize()
{
	for (int i = 1; i <= 8; i++)
		if (check(x, y, dx[i], dy[i]))
			arr[(x - 1) * m + y + dx[i] * m + dy[i]] = 1; //对于arr而言 1 1 其实是arr的第零个数字 arr共n * m 元素有效

	for (int i = 1; i <= n * m; i++)
		for (int j = 1; j <= n * m; j++)
			init.grid[i][j] = 0;
	for (int i = 1; i <= n * m; i++)
	{
		for (int k = 1; k <= 8; k++)
		{
			int ep1 = 0,ep2 = 0,ep3 = 1;
			if (i % m == 0) { ep1 = -1, ep2 = m, ep3 = 0; }
			if (check(i / m + ep3, i % m + ep2, dx[k], dy[k])) //在二维坐标下判断
				init.grid[i][m * (i / m + dx[k] + ep1)  + i % m + ep2 + dy[k]] = 1;        //这边每一行看作原来二维数组的一维形式
		}
	}
}

void output(int op)
{
	if (op == 1)
	{
		int cnt = 0;
		while (cnt < n * m)
		{
			cout << arr[cnt + 1] << " ";
			cnt++;
			if (cnt % m == 0)
				cout << endl;
		}
	}
	else if (op == 2)
	{
		int cnt = 0;
		while (cnt < n * m)
		{
			cout << ans[cnt + 1] << " ";
			cnt++;
			if (cnt % m == 0)
				cout << endl;
		}
	}
}

int main()
{
	cin >> n >> m >> x >> y >> k;

	initialize();

	if (k == 1) output(1);
	else
	{
		Mat a;
		a = qpow(init, k - 2);
		int cnt = 0;
		for (int i = 1; i <= n * m; i++)
			for (int j = 1; j <= n * m; j++)
				ans[i] = (ans[i] + arr[j] * a.grid[i][j]) % MOD;

		output(2);
	}
}