POJ - Matrix Power Series 【矩阵快速幂+二分求和】

Description

Given a n × n matrix A and a positive integer k, find the sum S = A + A² + A³ + … + A^k.

Input

The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 10⁹) and m (m < 10⁴). Then follow n lines each containing n nonnegative integers below 32,768, giving A’s elements in row-major order.

Output

Output the elements of S modulo m in the same way as A is given.

Sample Input

2 2 4
0 1
1 1

Sample Output

1 2
2 3

这题利用矩阵快速幂 + 二分求和，题目要求是Given a n × n matrix A and a positive integer k, find the sum S = A + A² + A³ + … + A^k.

⑴  这里，假设k = 6，可以得到 S = a + a^2 + a^3 + a^4 + a^5 + a^6 --->  S = (a + a^2 + a^3) + a^3(a^1 + a^2 + a^3) 

这里用A = (a + a^2 + a^3) ，那么，可以得到 S = A + a^3*A ，然后，然后，继续将A二分下去，直到A = a，利用递归求解就可以得到结果。

⑵ 这里还有一点要注意，当K的值为奇数或者偶数的时候，要分下面两种情况：

     ①  当k为偶数，假设k = 2n，S = a + a^2 + a^3+...+ a^n + (a + a^2 + a^3 +...+ a^n)a^n

     ②  当k为奇数，假设k = 2n+1，S = a + a^2 + a^3+...+ a^n + (a + a^2 + a^3 +...+ a^n)a^n + a^2n+1

下面是具体的AC代码：

#include <iostream>
#include <cstring>

using namespace std;
int size,mod,n; // size是n*n矩阵的大小，mod是取模大小，n是次数
class Matrix{
public:
	int data[30][30];
	Matrix();
	Matrix operator *(const Matrix rig);
	Matrix operator +(const Matrix rig);
	void showMatrix();  // 输出矩阵
	void unitMatrix();  // 单位矩阵
};
Matrix::Matrix()   // 矩阵的构造函数
{
	memset(data,0,sizeof(data));
}
Matrix Matrix::operator *(const Matrix rig)
{
	Matrix ans;
	for(int i = 0; i < size; i++)
	{
		for(int j = 0; j < size; j++)
		{
			for(int k = 0; k < size; k++)
			{
				ans.data[i][j] += data[i][k] * rig.data[k][j] % mod;
				ans.data[i][j] %= mod;
			}
		}
	}
	return ans;
}
Matrix Matrix::operator +(const Matrix rig)
{
	Matrix ans;
	for(int i = 0; i < size; i++)
	{
		for(int j = 0; j < size; j++)
		{
			ans.data[i][j] = (data[i][j] + rig.data[i][j]) % mod;
		}
	}
	return ans;
}
void Matrix::showMatrix()
{
	for(int i = 0; i < size; i++)
	{
		cout << data[i][0] % mod;
		for(int j = 1; j < size; j++)
		{
			cout << " " << data[i][j] % mod;
		}
		cout << '\n';
	}
}
void Matrix::unitMatrix()
{
	for(int i = 0; i < size; i++)
	{
		data[i][i] = 1;
	}
}

Matrix powMatrix(Matrix rig,int t)
{
	Matrix ans;
	ans.unitMatrix();
	while(t)
	{
		if(t & 1) ans = ans * rig;
		rig = rig * rig;
		t >>= 1;
	}
	return ans;
}
Matrix solve(Matrix rig, int t)
{
	Matrix tmp;
	tmp.unitMatrix();
	if(t == 1)
	{
		return rig;
	}
	else
	{
		if(t & 1)     // 相当于 A*a^t+A + a^(2t+1) ---> A*(a^t+1) + a^(2t+1)
			return (powMatrix(rig,t>>1)+ tmp)*solve(rig,t>>1) + powMatrix(rig,t);
		else          // 相当于 A*a^t+A ---> A*(a^t+1)
			return (powMatrix(rig,t>>1)+tmp)*solve(rig,t>>1);
	}
}
int main(int argc, char const *argv[])
{
	while(cin >> size >> n >> mod)
	{
		int a[size][size];
		Matrix res,s;
		for(int i = 0; i < size; i++)
		{
			for(int j = 0; j < size; j++)
			{
				cin >> s.data[i][j];
			}
		}
		res = solve(s,n);
		res.showMatrix();
	}
	return 0;
}