矩阵快速幂

矩阵快速幂

1. 什么是矩阵快速幂

在矩阵加速斐波那契数列中，我们往往需要求这个式子：

\[A^n \times \begin{bmatrix} f(1) \\ f(0) \end{bmatrix} = \begin{bmatrix} f(n+1) \\ f(n) \end{bmatrix} \]

如何快速的求一个矩阵\(A\)的\(k\)次幂呢？我们可以使用快速幂。即矩阵快速幂。
注：本文不涉及基本的线性代数知识。

2. 算法细节

我们使用结构体来方便的声明矩阵：

struct Matrix
{
	long long m[N][N];
}unit;

还需要重载*运算符来执行两个矩阵相乘的操作(Mod 为 1e9 + 7)：

Matrix operator * (const Matrix &x,const Matrix &y)
{
	Matrix ans;
	for(int i = 1;i <= n;i ++ )
	{
		for(int j = 1;j <= n;j ++ )
		{
			long long t = 0;
			for(int k = 1;k <= n;k ++ )
			{
				t += x.m[i][k] * y.m[k][j] % Mod;
			}
			ans.m[i][j] = t % Mod;
		}
	}
	return ans;
}

重载运算符之后套用快速幂模板得到矩阵快速幂的核心部分：

Matrix fastM(Matrix a,long long k)
{
	Matrix res = unit;
	
	while(k)
	{
		if(k & 1) res = res * a;
		a = a * a;
		k >>= 1;
	}
	return res;
}

你或许会发现上方的代码里有一句：

Matrix res = unit;

为什么呢？因为unit是一个单位矩阵，unit乘上任何矩阵都等于他本身。即\(n \times unit = n\)。如果我们不乘以单位矩阵，res将一直为0。
unit矩阵初始化：

void init()
{
	for(int i = 1; i <= n;i ++ )
	{
		unit.m[i][i] = 1;
	}
}

3. 完整代码

#include<bits/stdc++.h>
using namespace std;

const int N = 110;
const int Mod = 1e9 + 7;

long long n,k; 

struct Matrix
{
	long long m[N][N];
}unit;

Matrix operator * (const Matrix &x,const Matrix &y)
{
	Matrix ans;
	for(int i = 1;i <= n;i ++ )
	{
		for(int j = 1;j <= n;j ++ )
		{
			long long t = 0;
			for(int k = 1;k <= n;k ++ )
			{
				t += x.m[i][k] * y.m[k][j] % Mod;
			}
			ans.m[i][j] = t % Mod;
		}
	}
	return ans;
}

void init()
{
	for(int i = 1; i <= n;i ++ )
	{
		unit.m[i][i] = 1;
	}
}

Matrix fastM(Matrix a,long long k)
{
	Matrix res = unit;
	
	while(k)
	{
		if(k & 1) res = res * a;
		a = a * a;
		k >>= 1;
	}
	return res;
}

int main()
{
	Matrix a;
	cin >> n >> k;
	init();
	for(int i = 1;i <= n;i ++ )
	{
		for(int j = 1;j <= n;j ++ )
		{
			cin >> a.m[i][j];
		}
	}
	
	a = fastM(a,k);
	
	for(int i = 1;i <= n;i ++ )
	{
		for(int j = 1;j <= n;j ++ )
		{
			cout << a.m[i][j] << " ";
		}
		cout << endl;
	}
	return 0;
}

矩阵快速幂不难，在理解快速幂的基础上仔细思考就可理解。