矩阵快速幂

矩阵乘法在OI数论范畴内是一类很常见的的算法，用于优化线性递推方程。

定义

一个 \(m \times n\) 的矩阵 \(A\) 是一个由 \(m\) 行 \(n\) 列元素排列成的矩形阵列。即形如

\[A = \begin{bmatrix} a_{1 1} & a_{1 2} & \cdots & a_{1 n} \\ a_{2 1} & a_{2 2} & \cdots & a_{2 n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m 1} & a_{m 2} & \cdots & a_{m n} \end{bmatrix} \]

两个矩阵相乘的定义为一个 \(n\times t\) 的矩阵 \(A\) 和一个 \(t\times m\) 的矩阵 \(B\) 相乘，得到的结果为一个 \(n\times m\) 的矩阵 \(C\)。

其中 \(C\) 的形成遵循以下式子：

\[C_{i,j}=\sum_{k=1}^{t} (A_{i,k}\times B_{k,j}) \]

具体的，我们给出一个矩阵乘法的例子：

\[\left[\begin{matrix}1~~~~2~~~~3\\4~~~~5~~~~6\\\end{matrix}\right] \times \left[\begin{matrix}7 ~~~~~ 8\\9 ~~~~ 10\\11 ~~~ 12\\\end{matrix}\right]=\left[\begin{matrix}{1\times 7+2\times 9+3\times 11} ~~~~ {1\times 8+2\times 10 +3\times 12}\\{4\times 7+5\times 9+6\times 11} ~~~~ {4\times 8+5\times 10+6\times 12} \end{matrix}\right]=\left[\begin{matrix}58 ~~~~ 64\\139 ~~~~ 154\end{matrix}\right] \]

类似的，我们也可以定义出矩阵加法和矩阵减法，但这两种运算在竞赛中不常见，所以不讲了，再见。

两个矩阵若行列数不满足上述条件，则无法进行运算喵~

接下来给出一些性质和定义：

矩阵乘法遵循结合律和分配率，但不遵循交换律。

一个矩阵的乘方定义为多个该矩阵的乘积矩阵，根据定义，只有一个行列数相等的矩阵才能进行乘方运算。

特殊地，定义 \(n\times n\) 的单位元矩阵为：

\[E=\begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{bmatrix} \]

即对角线上的元素为 \(1\)，其余元素均为 \(0\)。

容易发现，任何一个列数为 \(n\) 的矩阵与 \(E\) 相乘，得到的结果还是这个矩阵，因此 \(E\) 的定义可以类比乘法中的 \(\times 1\) 来理解。

理解以后可以去做一做模板题来加深印象并形成代码模板。

这里作者推荐用结构体存储矩阵并对矩阵乘法进行封装，这样可以使代码更有可读性。

我的代码：

点击查看代码

struct no
{
	int num[50][50];//矩阵信息
	int n,m;	//矩阵行列数
	no(){memset(num,0,sizeof num);}//矩阵初始化
	int* operator [] (int i){return num[i]; }//方便书写
	inline friend no operator * (no a,no b)//封装的矩阵乘法
	{
		no c;
		c.n=a.n;c.m=b.m;
		for(int i=1;i<=c.n;i++)
			for(int j=1;j<=c.m;j++)
				for(int k=1;k<=a.m;k++)
					c[i][j]=(c[i][j]+a[i][k]*b[k][j]%mod)%mod;
		return c;
	}
};

应用

那么定义讲完了，这个东西究竟该怎么优化递推呢？

接下来就讲一讲矩阵乘法最普遍的应用————

矩阵快速幂

相信快速幂大家都学过，这里先放一下模板：

int qw(int x,int y)
{
	int ans=1ll;
	while(y)
	{
		if(y&1)ans=ans*x%mod;
		x=x*x%mod;
		y>>=1;
	}
	return ans*%mod;
}

那么现在对于两个矩阵 \(A\) 和 \(B\)，如果我们要计算 \(A^B\)，是否能用快速幂来加速呢？答案是肯定的！

如果用上文作者所说的方法存储矩阵，那么矩阵快速幂和一般快速幂的差别就只剩 \(ans\)，\(x\) 和函数的类型了（在进行乘法封装的时候已进行取模运算，所以在矩阵快速幂中无需取模）。

放一下模板（封装）：

inline friend no operator ^ (no x,int y)
{
	no ans;
	ans.n=ans.m=x.n;
	for(int i=1;i<=ans.n;i++)ans[i][i]=1;//初始为单位元矩阵
	while(y)
	{
		if(y&1) ans=ans*x;
		x=x*x;
		y>>=1;
	}
	return ans;
}

好了，现在基础都讲完了，那么这东西该怎么用呢？

题目

给你一个 \(n\)，让你输出 Fibonacci 数列第 \(n\) 项的值。答案对 \(998244353\) 取模。

Fibonacci 数列的递归定义如下：（就是斐波那契数列，作者为了高级用了英文）

\[Fib_i=\begin{cases}1~~~~ \text{if} ~~i=1 ~~\text{or} ~~i=2\\Fib_{i-1}+Fib_{i-2} ~~~~~~\text{else}\end{cases} \]

这是一道很简单的题，对吧？

但如果 \(n\leq 10^9\) 呢？

这时候我们就可以用矩阵快速幂来做这件事情啦！

怎么做呢？

首先我们给出 Fibonacci 数列的一般递推方程：

\[Fib_i=Fib_{i-1}+Fib_{i-2} \]

接下来我们做一件神奇的事情，把这个方程变成矩阵乘法的形式。

我们设一个状态矩阵为 \(\left[\begin{matrix}Fib_{i-1}~~~~Fib_{i-2}\\\end{matrix}\right]\)。

那么类似的，我们现在需要得到目标矩阵 \(\left[\begin{matrix}Fib_{i}~~~~Fib_{i-1}\\\end{matrix}\right]\)。

所以我们就需要构造一个转移矩阵，使得当前矩阵与它相乘以后可以得到目标状态的矩阵。

构造转移矩阵是一个矩阵乘法题目中最核心的思考部分。对于这道题，我们不难想到这个矩阵其实就是 \(\left[\begin{matrix}1~~~~1\\1~~~~0\\\end{matrix}\right]\)。

所以回到题目要求的东西，我们可以很自然地写出答案所在的矩阵 \(\left[\begin{matrix}Fib_{n+1}~~~~Fib_{n}\\\end{matrix}\right]\) 其实就是这个东西：

\[\left[\begin{matrix}Fib_2~~~~Fib_1\\\end{matrix}\right] \times \left[\begin{matrix}1~~~~1\\1~~~~0\\\end{matrix}\right]^{n-1} \]

那么我们可以神奇地发现，这玩意可以用矩阵快速幂来做！

所以我们可以愉快的给出代码了。

点击查看代码

#include<bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
#include<ext/pb_ds/hash_policy.hpp>
#include<ext/pb_ds/trie_policy.hpp>
#include<ext/pb_ds/priority_queue.hpp>
#define int long long
using namespace std;
using namespace  __gnu_pbds;
//gp_hash_table<string,int>mp2;
//__gnu_pbds::priority_queue<int,less<int>,pairing_heap_tag> q;
inline int read()
{
	int w=1,s=0;char ch=getchar();
	while(!isdigit(ch)){if(ch=='-')w=-1;ch=getchar();}
	while(isdigit(ch)){s=s*10+(ch-'0');ch=getchar();}
	return w*s;
}
const int mod=998244353;
const int maxn=1e6+10;
int k;
int b[maxn],c[maxn],sum[maxn];
int n,m; 
struct no
{
	int num[9][9];
	int n,m;	
	no(){memset(num,0,sizeof num);}
	int* operator [] (int i){return num[i];}
	inline friend no operator * (no a,no b)
	{
		no c;
		c.n=a.n;c.m=b.m;
		for(int i=1;i<=c.n;i++)
			for(int j=1;j<=c.m;j++)
				for(int k=1;k<=a.m;k++)
					c[i][j]=(c[i][j]+a[i][k]*b[k][j]%mod)%mod;
		return c;
	}
	inline friend no operator ^ (no x,int y)
	{
		no ans;
		ans.n=ans.m=x.n;
		for(int i=1;i<=ans.n;i++)
			ans[i][i]=1;
		while(y)
		{
			if(y&1) ans=ans*x;
			x=x*x;
			y>>=1;
		}
		return ans;
	}
};
signed main()
{
    cin>>n;
    no x,st;
    x.n=x.m=2;
    x[1][1]=x[1][2]=x[2][1]=1;//转移矩阵 
    st.n=1;st.m=2;
    st[1][1]=st[1][2]=1;//初始矩阵 
    no res;
    res=st*(x^(n-1));
    cout<<res[1][2];
    return 0;
}

本文仅介绍矩阵乘法，练习题请读者自行寻找喵（其实就是作者懒）。