【矩阵快速幂】封装类及测试用例及样例

作者推荐

视频算法专题

通俗的说，就是矩阵的乘方。

封装类

核心代码

class CMat
{
public:
	// 矩阵乘法
	static vector<vector> multiply(const vector<vector>& a, const vector<vector>& b) {
		const int r = a.size(), c = b.front().size(),iK = a.front().size();
		assert(iK == b.size());
		vector<vector> ret(r, vector(c));
		for (int i = 0; i < r; i++)
		{
			for (int j = 0; j < c ; j++) 
			{
				for (int k = 0; k < iK; k++)
				{
					ret[i][j] = (ret[i][j] + a[i][k] * b[k][j] ) % s_llMod;
				}
			}
		}
		return ret;
	}

	// 矩阵快速幂
	static vector<vector> pow( const vector<vector>& a, vector<vector> b, long long n) {
		vector<vector> res = a;
		for (; n; n /= 2) {
			if (n % 2) {
				res = multiply(res, b);
			}
			b = multiply(b, b);
		}
		return res;
	}
	static vector<vector> TotalRow(const vector<vector>& a)
	{
		vector<vector> b(a.front().size(), vector(1, 1));
		return multiply(a, b);
	}
protected:
	const static long long s_llMod = 1e9 + 7;
};

测试用例

template
void Assert(const T& t1, const T& t2)
{
	assert(t1 == t2);
}

template
void Assert(const vector& v1, const vector& v2)
{
	if (v1.size() != v2.size())
	{
		assert(false);
		return;
	}
	for (int i = 0; i < v1.size(); i++)
	{
		Assert(v1[i], v2[i]);
	}
}


int main()
{
	vector<vector> pre = { {1,2} };
	vector<vector> mat = { {2,3},{1,10} };
	{	
		auto res = CMat::pow(pre, mat, 0);
		Assert(pre, res);
	}
	{
		auto res = CMat::multiply(pre, mat);
		Assert(vector<vector>{ {4, 23}}, res);
		auto res2 = CMat::pow(pre, mat,1);
		Assert(res2, res);
	}
	{
		auto res = CMat::pow(pre, mat, 2);
		auto res1 = CMat::multiply(pre, mat);
		auto res2 = CMat::multiply(res1, mat);
		Assert(res2, res);
		Assert(vector<vector>{ {31, 242}}, res);
	};

	for (int i = 3; i < 100; i++)
	{
		auto res = pre;
		for (int j = 0; j < i; j++)
		{
			res = CMat::multiply(res, mat);
		}
		auto res2 = CMat::pow(pre, mat, i);
		Assert(res2, res);
	}
	
}

具体例子

题目、分析和原理见：

【动态规划】【矩阵快速幂】【滚动向量】C++算法552. 学生出勤记录 II

原解法用二维表示状态，改成一维。 i是缺勤数量,j是连续迟到数，新的状态为：3*i+j

6种状态，故转移矩阵为6行6列，故结果矩阵为6列，6个数据1行就足够了。

令旧结果矩阵为mat1，转移矩阵为mat2，新矩阵为mat3，K mat1的列数，mat2的行数。则：

mat3[r][c] = Sum

[

)

^{i}_{[0,k)}

[0,k)i(mat1[r][i]*mat2[i][c])

i在mat1中列号，在mat2中是行号。也就是旧状态在第几列，mat2就在第几行。

新状态就是mat2的行号。

class Solution {
public:
	int checkRecord(int n) {
		vector<vector> pre(1, vector(6));//1行6列	
		pre[0][0] = 1;
		vector<vector> mat(6, vector(6));
		{	
			//之前的状态在pre是第几列，矩阵中就是第几行。新状态的列号就矩阵的列号
			//处理一次缺勤 ，缺勤两次排除
			for (int i = 0; i < 3; i++)
			{
				mat[i][3]++;
			}
			//处理请假
			for (int i = 0; i < 2; i++)
			{
				for (int j = 0; j < 2; j++)
				{
					const int pre = 3 * i + j;
					mat[pre][pre + 1]++;
				}
			}
			//处理正常
			for (int i = 0; i < 2; i++)
			{
				for (int j = 0; j < 3; j++)
				{
					const int pre = 3 * i + j;
					const int cur = 3 * i;
					mat[pre][cur]++;
				}
			}
		}
		auto res = CMat::pow(pre, mat, n);
		res = CMat::TotalRow(res);
		return res[0][0];
	}
};

测试用例

template

void Assert(const T& t1, const T& t2)

{

assert(t1 == t2);

}

template

void Assert(const vector& v1, const vector& v2)

{

if (v1.size() != v2.size())

{

assert(false);

return;

}

for (int i = 0; i < v1.size(); i++)

{

Assert(v1[i], v2[i]);

}

int main()

{

int n;

{

Solution sln;

n = 0;