矩阵快速幂 超详细介绍

矩阵概念

什么是矩阵

在数学中，矩阵（Matrix）是一个按照长方阵列排列的复数或实数集合

特殊矩阵（单位矩阵）

单位矩阵之所以是特殊矩阵是因为单位矩阵的主对角线上全是1，其他元素全为0，用单位矩阵乘以任何矩阵都将得到原矩阵，类似数字1的效果。
下面是一个3x3的单位矩阵

矩阵运算

矩阵加法减法比较简单这里就不过多介绍了，主要说一下矩阵乘法，

矩阵乘法

矩阵相乘只有在第一个矩阵的列数（column）和第二个矩阵的行数（row）相同时才有意义，设A为m x p 的矩阵，B为p x n 的矩阵，那么称m x n的矩阵C为矩阵A与B的乘积 ，其中矩阵C中的第 行第 列元素可以表示为

矩阵是满组结合律的，而就是因为矩阵满足了结合律，所以矩阵乘法可以使用快速幂算法，快速幂算法中我们把一个数连续乘的问题利用结合律简化，使运算循环次数减少，例如：计算3⁵,原本要计算的是3 x 3 x 3 x 3 x 3，但是用结合律简化得（3 x 3）x （3 x 3）x 3进而得9 x 9 x 3 继续使用结合律（9 x 9）x 3 最后得81 x 3，把需要五次累乘问题转换成了一次累乘问题，在幂次小的时候没多少区别，但是幂次大的时候差距会非常明显后者更有效率，关于快速幂思想就不再过多说了，对此不明白的同学可以参考参考专门写的一篇快速幂的博客，那里写的很详细点我传送。矩阵因为有结合律也可以用次方法简化，例如M是一个矩阵现在要求M⁷。传统方法是七次连续乘M也就是MxMxMxMxMxMxM 。但根据快速幂思想加上矩阵支持结合律可以这样写（MxM）x（MxM）x（MxM）xM也就是求M²的3次方然后再乘以一个M得M² x M² x M² x M我们可以继续用结合律（M² x M² ）x M² x M得M⁴x M² x M,这就是快速幂思想在矩阵上的运用。矩阵快速幂的实现代码其实就是基于快速幂代码模板的。
矩阵快速幂样例代码

矩阵快速幂代码C++实现

#include<iostream>

using namespace std;

const int MAX = 15;
const int MOD = 9973;

struct Matrix
{
	int mat[MAX][MAX];
};

Matrix create_identity_matrix(Matrix E, int n);		//构造单位矩阵
Matrix mat_mul(Matrix a, Matrix b, int n);			//矩阵相乘
Matrix quick_power(Matrix basic, int power, int n);	//矩阵快速幂
void input_mat(Matrix* p, int n, int m);			//输入矩阵
void print_marrix(Matrix a, int n, int m);			//打印矩阵

int main(void)
{
	Matrix mat_now;
	Matrix ans_now;
	int n, k;
	//n是矩阵的行列数， k是幂次数
	cin >> n >> k;
	input_mat(&mat_now, n, n);
	ans_now = quick_power(mat_now, k, n);
	print_marrix(ans_now, n, n);

	return 0;
}

Matrix create_identity_matrix(Matrix E, int n)
{
	for (int i = 0; i < n; i++)
		for (int j = 0; j < n; j++)
			if (i == j)
				E.mat[i][j] = 1;
			else
				E.mat[i][j] = 0;
	return E;
}

Matrix mat_mul(Matrix a, Matrix b, int n)
{
	Matrix ans;
	memset(ans.mat, 0, sizeof(ans.mat));
	//这里会看出三次循环顺序是i，k，j而不是说i，j，k，当然i， j， k也对
	//不过矩阵在计算机内存里是线性存储的所有i，k，j更快
	//也就是说i，j，k顺序计算适合人类计算而i，k，j更适合计算机计算
	for (int i = 0; i < n; i++)
		for (int k = 0; k < n; k++)
			for (int j = 0; j < n; j++)
			{
				ans.mat[i][j] += a.mat[i][k] * b.mat[k][j];
				ans.mat[i][j] %= MOD;
				//ans.mat[i][j] += a.mat[i][k] * b.mat[k][j] % MOD;
			}
	return ans;
}


Matrix quick_power(Matrix basic, int power, int n)
{
	Matrix result;
	//下面这一行本质上没有意义
	//但是vs2019会以为你调用了未初始化变量然后报错
	//当然这个机制可以在设置关闭，如果没有关闭的话加上下面这一行不然无法通过编译
	memset(result.mat, 0, sizeof(result.mat));
	result = create_identity_matrix(result, n);	//初始化为单位矩阵
	while (power > 0)
	{
		if (power & 1)
			result = mat_mul(result, basic, n);
		basic = mat_mul(basic, basic, n);
		power >>= 1;
	}
	return result;
}

void input_mat(Matrix* p, int n, int m)
{
	for (int i = 0; i < n; i++)
		for (int j = 0; j < m; j++)
			cin >> p->mat[i][j];
}

void print_marrix(Matrix a, int n, int m)
{
	for (int i = 0; i < n; i++)
	{
		for (int j = 0; j < m; j++)
			cout << ' ' << a.mat[i][j];
		cout << endl;
	}
}

代码分解

输入输出矩阵的函数

这个函数利用两个for循环输入与打印矩阵的。

void input_mat(Matrix* p, int n, int m)		//输入函数
{
	for (int i = 0; i < n; i++)
		for (int j = 0; j < m; j++)
			cin >> p->mat[i][j];
}

void print_marrix(Matrix a, int n, int m)	//输出函数
{
	for (int i = 0; i < n; i++)
	{
		for (int j = 0; j < m; j++)
			cout << ' ' << a.mat[i][j];
		cout << endl;
	}
}

矩阵相乘

因为矩阵乘法无法用内部的 * 运算符计算乘法，所有我们用函数实现相乘，注意里面实现取余操作的位置与方法。

Matrix mat_mul(Matrix a, Matrix b, int n)
{
	Matrix ans;
	memset(ans.mat, 0, sizeof(ans.mat));
	//这里会看出三次循环顺序是i，k，j而不是说i，j，k，当然i， j， k也对
	//不过矩阵在计算机内存里是线性存储的所有i，k，j更快
	//也就是说i，j，k顺序计算适合人类计算而i，k，j更适合计算机计算
	for (int i = 0; i < n; i++)
		for (int k = 0; k < n; k++)
			for (int j = 0; j < n; j++)
			{
				ans.mat[i][j] += a.mat[i][k] * b.mat[k][j];
				ans.mat[i][j] %= MOD;
				//快速幂中常常带有取余某个数的要求，在矩阵快速幂中这个取余需要放在乘法函数里
				//ans.mat[i][j] += a.mat[i][k] * b.mat[k][j] % MOD;
			}
	return ans;
}

构造单位矩阵函数

在快速幂中我们需要把存储答案的变量初始化为1，但是矩阵快速幂中不能是1因为矩阵不可以。代替1的是单位矩阵，我们要把存储答案的变量初始化为单位矩阵。

Matrix create_identity_matrix(Matrix E, int n)
{
	for (int i = 0; i < n; i++)
		for (int j = 0; j < n; j++)
			if (i == j)
				E.mat[i][j] = 1;
			else
				E.mat[i][j] = 0;
	return E;
}

矩阵快速幂实现函数

矩阵快速幂函数实现模板其实就是快速幂模板可以和快速幂模板代码做对比查看理解

Matrix quick_power(Matrix basic, int power, int n)
{
	Matrix result;
	//下面这一行本质上没有意义
	//但是vs2019会以为你调用了未初始化变量然后报错
	//当然这个机制可以在设置关闭，如果没有关闭的话加上下面这一行不然无法通过编译
	memset(result.mat, 0, sizeof(result.mat));
	result = create_identity_matrix(result, n);	//初始化为单位矩阵
	while (power > 0)
	{
		if (power & 1)
			result = mat_mul(result, basic, n);
		basic = mat_mul(basic, basic, n);
		power >>= 1;
	}
	return result;
}

快速幂模板代码

long long int quik_power(int base, int power, int p)
{
	long long int result = 1;   
	while (power > 0)           
	{
		if (power & 1)         							
			result = result * base % p;   
		base = base * base % p ;    						
		power >>= 1;         						
	}
	return result % p;       
}

矩阵快速幂应用

经典题目

题目HDOJ 1575 Tr A
题解 HDOJ 1575 Tr A 题解