首页 > 代码库 > poj 3233 Matrix Power Series

poj 3233 Matrix Power Series

Matrix Power Series
Time Limit: 3000MS Memory Limit: 131072K
Total Submissions: 15997 Accepted: 6840

Description

Given a n × n matrix A and a positive integer k, find the sum S = A + A2 + A3 + … + Ak.

Input

The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 109) and m (m < 104). 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

题意:给出矩阵A,求S = A + A^2 + A^3 + … + A^k ;

解题思路:二分+快速矩阵幂

在这里所谓的二分的意思是:举例说明假设k=18,则

S=A+A^2+A^3+A^4+A^5+A^6+A^7+A^8+A^9+A^10+A^11+A^12+A^13+A^14+A^15+A^16+A^17+A^18

S=(A+A^2+A^3+A^4+A^5+A^6+A^7+A^8+A^9)*(1+A^9)

S=(A+A^2+A^3+A^4)*(1+A^4)*(1+A^9)

S=(A+A^2)*(1+A^2)*(1+A^4)*(1+A^9)

S=A*(1+A)*(1+A^2)*(1+A^4)*(1+A^9)

参考代码如下:

#include <iostream>
#include <string.h>
using namespace std;
int n,m,k;
struct Matrix{
	int mat[31][31];
};
Matrix add(Matrix a,Matrix b){
	Matrix ans;
	for (int i=0;i<n;i++){
		for (int j=0;j<n;j++){
			ans.mat[i][j]=a.mat[i][j]+b.mat[i][j];
			if (ans.mat[i][j]>=m){
				ans.mat[i][j]%=m;
			}
		}
	}
	return ans;
}
Matrix mul(Matrix a,Matrix b){
	Matrix ans;
	for (int i=0;i<n;i++){
		for (int j=0;j<n;j++){
			ans.mat[i][j]=0;
			for (int k=0;k<n;k++){
				ans.mat[i][j]+=a.mat[i][k]*b.mat[k][j];
				if (ans.mat[i][j]>=m){
					ans.mat[i][j]%=m;
				}
			}
		}
	}
	return ans;
}
Matrix Init(){
	Matrix ans;
	for (int i=0;i<n;i++){
		for (int j=0;j<n;j++){
			if (i==j)
				ans.mat[i][j]=1;
			else
				ans.mat[i][j]=0;
		}
	}
	return ans;
}
Matrix exp(Matrix a,int k){
	Matrix ans=Init();
	while (k){
		if (k&1)
			ans=mul(ans,a);
		a=mul(a,a);
		k>>=1;
	}
	return ans;
}
Matrix Calculate(Matrix a,int k){
	Matrix ans;
	memset(ans.mat,0,sizeof(ans.mat));
	Matrix temp=Init();
	if (k==1)
		return a;
	if (k%2)
		ans=add(Calculate(a,k-1),exp(a,k));
	else{
		Matrix s=Calculate(a,k/2);
		ans=add(s,mul(s,exp(a,k/2)));
	}
	return ans;
}
int main(){
	Matrix a;
	while (cin>>n>>k>>m){
		for (int i=0;i<n;i++){
			for (int j=0;j<n;j++){
				cin>>a.mat[i][j];
				if (a.mat[i][j]>=m){
					a.mat[i][j]%=m;
				}
			}
		}
		Matrix ans=Calculate(a,k);
		for (int i=0;i<n;i++){
			for (int j=0;j<n-1;j++)
				cout<<ans.mat[i][j]<<" ";
			cout<<ans.mat[i][n-1]<<endl;
		}
	}
	return 0;
}





poj 3233 Matrix Power Series