矩阵快速幂运用

引入

想要了解矩阵快速幂，就不得不提到矩阵的概念。矩阵就像一个二维数组，存储了一组数据，如：

$M=\left[ \begin{matrix} 1&2&3\\ 4&5&6\\ 7&8&9\end{matrix} \right]$

用 $M_{i,j}$ 表示矩阵第 $i$ 行第 $j$ 列的数据，如 $M_{2,3}=6$。

矩阵运算

加法

矩阵的加法和实数加法类似，要求运算的两个矩阵大小相等，将对应位置相加即可。

减法

与加法相同，对应位置相减即可。

乘法

矩阵乘法并非对应位置相乘。

若矩阵 $A\times B=C$ 有意义，那么 $A$ 的列数要求与 $B$ 的行数相等。

若 $A$ 的列数与 $B$ 的行数等于 $n$：

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

即结果第 $i$ 行第 $j$ 列的值为 $A$ 的第 $i$ 行与 $B$ 的第 $j$ 列的对应项相乘后求和。

一个很直观的理解是，把矩阵 $A$ 放在结果的左侧， 矩阵 $B$ 放在结果的下方，将行与列延长，对应相乘后结果写在交点处。

$\left[ \begin{matrix} \bold{1}&\bold{2}&\bold{3}\\ 4&5&6\\ 7&8&9\end{matrix} \right]\left[ \begin{matrix} \rightarrow&36&&\\ &\uparrow&&\\ &\uparrow&&\end{matrix} \right]\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \left[ \begin{matrix} 1&\bold{2}&3\\ 4&\bold{5}&6\\ 7&\bold{8}&9\end{matrix} \right]$

注意矩阵乘法满足结合律，但是不一定满足交换律，$A\times B$ 成立不一定代表 $B\times A$ 成立。

除法

$\rightarrow$To be continued.

快速幂

快速幂可以 $O(\log n)$ 求解 $a^n$，对矩阵同样适用。
代码：

int ans=1;
while(n){
	if(n&1)ans*=a;
	a*=a;
	n>>=1;
}

矩阵快速幂

可以使用重载运算符的办法直接进行矩阵快速幂，方法和普通快速幂一样，但是要注意乘法顺序。ans的初始值要特别注意。如果想让它充当1，则需要设置成单位矩阵 $I$。

$\left[\begin{matrix}1&0&0&\cdots&0\\0&1&0&\cdots&0\\0&0&1&\cdots&0\\\vdots&\vdots&\vdots&\ddots&0\\0&0&0&0&1\end{matrix}\right]$

递推

使用矩阵快速幂，可以用 $O(\log n)$ 的时间复杂度完成 $n$ 次递推关系，如斐波那契数列。矩阵可以通过递推关系推算出来。

$\because\left[\begin{matrix}1&1\\0&1\end{matrix}\right] \times \left[\begin{matrix}F_{n-1}\\F_{n-2}\end{matrix}\right] = \left[\begin{matrix}F_n\\F_{n-1}\end{matrix}\right]\\$

$\therefore F_n=\left[\begin{matrix}1&1\\0&1\end{matrix}\right]^n \times \left[\begin{matrix}1\\1\end{matrix}\right]$

推算矩阵

矩阵由转移关系得来。以斐波那契数列为例：

$\because F_n=1F_{n-1}+1F_{n-2}$

$F_{n-1}=1F_{n-1}+0F_{n-2}$

$\therefore M=\left[\begin{matrix}1&1\\0&1\end{matrix}\right]$

其余类似。

例题

$F(n,1)=F(n-1,2)+2F(n-3,1)+5$

$F(n,2)=F(n-1,1)+3F(n-3,1)+2F(n-3,2)+3$

已知

$F(1,1)=2,F(1,2)=3,F(2,1)=1,F(2,2)=4,F(3,1)=6,F(3,2)=5$

求 $F(n,1)$ 和 $F(n,2)$。

矩阵

$\left[ \begin{matrix} 0&1&0&0&2&0&5\\ 1&0&0&0&3&2&3\\ 1&0&0&0&0&0&0\\ 0&1&0&0&0&0&0\\ 0&0&1&0&0&0&0\\ 0&0&0&1&0&0&0\\ 0&0&0&0&0&0&1\\ \end{matrix} \right]$

代码

#include "bits/stdc++.h"
using namespace std;
const long long mod=99999999;
long long base[7][7]={
	{0,1,0,0,2,0,5},
	{1,0,0,0,3,2,3},
	{1,0,0,0,0,0,0},
	{0,1,0,0,0,0,0},
	{0,0,1,0,0,0,0},
	{0,0,0,1,0,0,0},
	{0,0,0,0,0,0,1},
};
struct matrix{
	long long v[30][30];
	long long x,y;
};
matrix operator+(matrix a,matrix b){
	matrix c;
	if(a.x!=b.x||a.y!=b.y)return c;
	long long x=a.x,y=a.y;
	c.x=x;c.y=y;
	for(long long i=0;i<x;i++){
		for(long long j=0;j<y;j++){
			c.v[i][j]=a.v[i][j]+b.v[i][j];
		}
	}
	return c;
}
matrix operator-(matrix a,matrix b){
	matrix c;
	if(a.x!=b.x||a.y!=b.y)return c;
	long long x=a.x,y=a.y;
	c.x=x;c.y=y;
	for(long long i=0;i<x;i++){
		for(long long j=0;j<y;j++){
			c.v[i][j]=a.v[i][j]-b.v[i][j];
		}
	}
	return c;
}
matrix operator*(matrix a,matrix b){
	matrix c;
	if(a.y!=b.x)return c;
	long long x=a.x,y=b.y,z=a.y;
	c.x=x;c.y=y;
	for(long long i=0;i<x;i++){
		for(long long j=0;j<y;j++){
			c.v[i][j]=0;
			for(long long k=0;k<z;k++){
				c.v[i][j]=(c.v[i][j]+a.v[i][k]*b.v[k][j]+mod)%mod;
			}
		}
	}
	return c;
}
void init(matrix &a){
	for(long long i=0;i<a.x;i++){
		for(long long j=0;j<a.y;j++){
			a.v[i][j]=base[i][j];
		}
	}
}
ostream& operator<<(ostream& ous,matrix a){
	for(long long i=0;i<a.x;i++){
		for(long long j=0;j<a.y;j++){
			ous<<a.v[i][j]<<' ';
		}
		ous<<endl;
	}
	return ous;
}
long long n;
matrix& operator *=(matrix &a,matrix n){
	a=a*n;
	return a;
}
matrix a,ans;
int main(){
	cin>>n;
	a.x=a.y=7;
	init(a);
	ans.x=7;ans.y=1;
	ans.v[0][0]=6;
	ans.v[1][0]=5;
	ans.v[2][0]=1;
	ans.v[3][0]=4;
	ans.v[4][0]=2;
	ans.v[5][0]=3;
	ans.v[6][0]=1;
	if(n==1){
		cout<<2<<endl<<3;
	}else if(n==2){
		cout<<1<<endl<<4;
	}else if(n==3){
		cout<<6<<endl<<5;
	}
	if(n<4)return 0;
	n-=3;
	while(n){
		if(n&1)ans=a*ans;
		a*=a;
		n>>=1;
	}
	cout<<ans.v[0][0]<<endl<<ans.v[1][0];
	return 0;
}