poj3233——Matrix Power Series

矩阵乘法。

S = A + A² + A³ + … + A^k ，把问题转化以加速，令

B =A I

0 I

则B^(k + 1) = A^(k + 1) I + A + A² + A³ + … + A^k

0I

用二分法求B^(k + 1)

来源：http://archive.cnblogs.com/a/1960189/

#include<iostream>
#include<cstdio>
#include<cstring>
using namespace std;
class  node 
{
public:
	int matrix[66][66];
}ans,b;
int n,m,k;
void init(int n,int m)
{
	int i,j;
	memset(ans.matrix ,0,sizeof(ans.matrix ));
	memset(b.matrix ,0,sizeof(b.matrix ));
	for(i=1;i<=n;i++)
		for(j=1;j<=n;j++)
		{
			scanf("%d",&ans.matrix [i][j]);
			ans.matrix [i][j]=ans.matrix [i][j]%m;
			b.matrix [i][j]=ans.matrix [i][j];
		}
		for(i=1;i<=n;i++)  {ans.matrix [i][n+i]=1; b.matrix [i][n+i]=1;}
		for(i=1+n;i<=n*2;i++) { ans.matrix [i][i]=1;b.matrix [i][i]=1;} 
}
void mul(node &a,node &b)
{
	node x;
	int i,k,j,sum;
	for(i=1;i<=n*2;i++)
		for(j=1;j<=n*2;j++)
			x.matrix[i][j]=a.matrix[i][j];
	for(i=1;i<=n*2;i++)
		for(j=1;j<=n*2;j++)
		{
			sum=0;
			for(k=1;k<=n*2;k++)
				sum+=(x.matrix [i][k]*b.matrix [k][j]);
			a.matrix [i][j]=sum%m;
		}
}
void solve()
{
	node a;
	int i,j;
	for(i=1;i<=n*2;i++)
		for(j=1;j<=n*2;j++)
		{
			a.matrix [i][j]=ans.matrix [i][j];
		}
		mul(ans,a);
}
void work(int num)
{
	if(num==1) return ;
	work(num/2);
	solve();
	if(num%2==1)
		mul(ans,b);
}
int main()
{
	int i,j;
	cin>>n>>k>>m;
	init(n,m);
	work(k+1);
	for(i=1;i<=n;i++)
	{
		if(ans.matrix [i][i+n]==0)
			ans.matrix [i][i+n]=m;
		ans.matrix [i][i+n]-=1;
	}
	for(i=1;i<=n;i++)
	{
		cout<<ans.matrix [i][n+1];
		for(j=2+n;j<=n*2;j++)
			cout<<" "<<ans.matrix [i][j];
		cout<<endl;
	}
	return 0;
}