poj 3233 Matrix Power Series

A为一个n*n的矩阵，求A+A^2+A^3+...+A^n

 

S_k = A + A² + A³ + … + A^k  

    =(1+A^k/2)*(A + A² + A³ + … + A^k/2  )+{A^k}

    =(1+A^k/2)*(S_k/2 )+{A^k}// k为偶数时无 {A^k}

A^k   可用二分迭代求出

因此,只要求出 上面的三部分就可以求出 S_k

设f(n)=A+A^2+A^3+...+A^n

n%2==1时，f(n)=f(n-1)+A^n

n%2==0时，f(n)=f(n/2)+f(n/2)*A^(n/2)

 

由于矩阵乘法满足结合律，计算A^n时，也可以二分

#include<iostream>
#include<cstdio>
using namespace std;
int n,m;
struct Mat
{
	int mat[31][31];
	Mat operator*(const Mat &x)
	{
		Mat tmp;
		int i,j,k;
		for(i=0;i<n;i++)
		{
			for(j=0;j<n;j++)
			{
				tmp.mat[i][j]=0;
				for(k=0;k<n;k++)
				{
					tmp.mat[i][j]+=(mat[i][k]*x.mat[k][j])%m;
					tmp.mat[i][j]%=m;
				}

			}
		}
		return tmp;
	}
	Mat operator+(const Mat &x)
	{
		Mat tmp;
		int i,j;
		for(i=0;i<n;i++)
			for(j=0;j<n;j++)
			{
				tmp.mat[i][j]=(mat[i][j]+x.mat[i][j])%m;
			}
		return tmp;
	}
}p;
Mat _pow(int k)
{
	if(k==1)
		return p;
	if(k&1)
		return _pow(k-1)*p;
	else
	{
		Mat tmp=_pow(k/2);
	
		return tmp*tmp;
	}
}
Mat cal(int k)
{
	if(k==1)
		return p;
	else
	{
		if(k&1)
			return cal(k-1)+_pow(k);
		else
		{
			Mat tmp=cal(k/2);
			return tmp+tmp*_pow(k/2);
		}
	}
}
int main()
{
	int i,j,k;
	scanf("%d%d%d",&n,&k,&m);
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
			scanf("%d",&p.mat[i][j]);
	}
	Mat res=cal(k);
	for(i=0;i<n;i++)
	{
		for(j=0;j<n-1;j++)
		{
			printf("%d ",res.mat[i][j]);
		}
		printf("%d\n",res.mat[i][j]);
	}
	return 0;
}