POJ 3233 Matrix Power Series(求矩阵幂的和——分块矩阵快速幂 or 二分递归+矩阵快速幂)

Matrix Power Series

Time Limit: 3000MS		Memory Limit: 131072K
Total Submissions: 21451		Accepted: 8993

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

题意：给定 n*n 的矩阵A和正整数 k 和 m 。 求矩阵A的幂之和。 S = A + A^2 + A^3 + ... + A^k    输出S的各元素对m取余的结果。  

知识储备： 矩阵   矩阵乘法  单位矩阵  分块矩阵  

题解：

方法一：构造法

因为分块矩阵依然满足矩阵乘法，所以可以构造成分块矩阵幂求解，详见 《挑战》 P.205

代码如下：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <vector>
#include <algorithm>
using namespace std;
#define LL long long 
#define maxn 35
typedef vector<LL> vec;
typedef vector<vec> mat;
int m;

mat mul(mat &A, mat &B)//矩阵乘法 
{
	mat C(A.size(), vec(B[0].size()));
	for(int i=0; i<A.size(); ++i)
	{
		for(int k=0; k<B.size(); ++k)
		{
			for(int j=0; j<B[0].size(); ++j)
				C[i][j] = (C[i][j] + A[i][k]*B[k][j]) % m;
		}
	}
	return C;
}

mat quick_pow(mat A,int n)//快速幂 
{
	mat B(A.size(), vec(A.size()));
	for(int i=0; i<A.size(); ++i)
		B[i][i]=1;
	while(n>0)
	{
		if(n&1)
			B = mul(B, A);
		A = mul(A, A);
		n >>= 1;
	}
	return B;
}

int main()
{
	int n,k;
	LL a;
	while(scanf("%d%d%d",&n,&k,&m)!=EOF)
	{
		mat A(n*2, vec(n*2));
		for(int i=0; i<n; ++i)
		{
			for(int j=0; j<n; ++j)
			{
				scanf("%lld",&a);
				A[i][j] = a;
			}
			A[n+i][i] = 1;	A[n+i][n+i] = 1;//构造出分块矩阵下方的两个单位矩阵块 
		}
		A = quick_pow(A, k+1);// I+A+A^2+...+A^k，构造后的 S(k) = I+A+...+A^k-1 ，所以这里我们要求出S(k+1) 
		for(int i=0; i<n; ++i)
		{
			for(int j=0; j<n ;++j)
			{
				a = A[n+i][j] % m;
				if(i==j)//减去最早加上的单位矩阵 
					a = (a+m-1) % m;
				if(j==n-1)
					printf("%lld\n",a);
				else
					printf("%lld ",a);
			}	
		}
	}
	return 0;
} 

方法二： 二分+矩阵快速幂

我们知道对于矩阵 A^x 可以用矩阵快速幂 快速的求解。 我们要考虑的是怎么减少 求解 S = A + A^2 + A^3 +...+A^k 过程中的 相加 这个操作的耗时。 下面我们把 矩阵A 看成一个普通的 数字A 

对于S(10)我们有： 

S(10) = A + A^2 + A^3 +A^4 + A^5 + A^5 * (A +  A^2 +A^3 +A^4 + A^5 )

对于S(5)我们有：

S(5) = A + A^2 +A^3 + A^3 * (A + A^2)

对于S(2)我们有：

S(2) = A + A * (A)

题中已经给出了 A 的值， 那么我们可以通过递归回溯得到S(2)的值，求出A^3的值，就可以通过上面的式子得到S(5)的值， 进而得到S(10)的值。 那么对于每个S(k)我们都可以这样二分递归求解。 注意这里的S(n)当n为奇数需要多加上一次A^(n/2+1)的值（ 例如上面的S(5) ）。 最后加上矩阵的快速幂 与 加法操作就行了。

代码如下：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
#define LL long long
#define maxn 40
int n,m,k;
struct Node
{
	LL mat[maxn][maxn];
}A;

Node mul(Node B, Node C)//矩阵乘法 
{
	Node D;
	for(int i=0; i<n; ++i)
		for(int j=0; j<n; ++j)
			D.mat[i][j] = 0;
	for(int i=0; i<n ;++i)
	{
		for(int k=0; k<n; ++k)
		{
			for(int j=0; j<n; ++j)
				D.mat[i][j] = (D.mat[i][j] + B.mat[i][k]*C.mat[k][j]) % m;
		}
	}
	return D;
}

Node quick_pow(Node B,int x)//快速幂 
{
	Node C;
	for(int i=0; i<n; ++i)
		for(int j=0; j<n; ++j)//初始化为单位矩阵 
			C.mat[i][j] = (i==j?1:0);
	while(x>0)
	{
		if(x&1)
			C = mul(C, B);
		B = mul(B, B);
		x >>= 1;
	}
	return C;
}

Node Add(Node B, Node C)//矩阵加法 
{
	Node D;
	for(int i=0; i<n; ++i)
	{
		for(int j=0; j<n; ++j)
			D.mat[i][j] = (B.mat[i][j] + C.mat[i][j]) % m;
	}
	return D;
}


Node dfs(int x)
{
	if(x==1)	
		return A;
	Node ans = dfs(x/2), temp;
	if(x&1){
		temp = quick_pow(A, x/2+1);
		ans = Add(ans, mul(temp, ans));
		ans = Add(temp, ans);//奇数情况要加上中间的矩阵 
	}
	else{
		temp = quick_pow(A, x/2);
		ans = Add(ans, mul(temp, ans));
	}
	return ans;
}

int main()
{
	while(scanf("%d%d%d",&n,&k,&m)!=EOF)
	{
		for(int i=0; i<n; ++i)
		{
			for(int j=0; j<n; ++j)
				scanf("%lld",&A.mat[i][j]);
		}
		Node ans = dfs(k);
		for(int i=0; i<n; ++i)
		{
			for(int j=0; j<n-1; ++j)
				printf("%lld ",ans.mat[i][j]);
			printf("%lld\n",ans.mat[i][n-1]);
		}	
	}
	return 0;
}