poj 3233 Matrix Power Series

分类：分治 难度：1.5

 

题意：给定矩阵a，求s=a^1+…+a^k，结果mod m。矩阵乘法，两次分治，计算a^k和a^1+…+a^k。即a^k=(a*a)^(k/2), a^1+...+a^k=a^1+...a^(k/2) + a^(k/2) * [ a^1+...a^(k/2) ]。

#include<cstdio>
#include<cstring>

int n,k,m;
const int N = 32;


typedef struct 
{
	int d[N][N];
}ma;
ma e;

ma madd(ma x, ma y)
{
	int i,j;
	ma z;
	for(i=0; i<n; i++)
	{
		for(j=0; j<n; j++)
		{
			z.d[i][j] = x.d[i][j] + y.d[i][j];
			z.d[i][j] %= m;
		}		
	}
	return z;
}

ma mmult(ma x, ma y)
{
	int i,j,l;
	int ans[N][N];
	ma z;
	for(i=0; i<n; i++)
	{
		for(j=0; j<n; j++)
		{
			z.d[i][j] = 0;
			for(l=0; l<n; l++)
			{
				z.d[i][j] += x.d[i][l]*y.d[l][j];
				z.d[i][j] %= m;
			}
		}
	}	
	return z;
}

ma mpow(ma a,int k)
{
	if(k==0)
		return e;
	if(k==1)
		return a;
	if(k%2)
		return mmult(mpow(a,k-1),a);
	return mpow(mmult(a,a),k/2);
}

void print(ma x)
{
	for(int i=0; i<n; i++)
	{
		for(int j=0; j<n; j++)
			printf("%d ",x.d[i][j]);
		printf("\n");
	}
}

ma doit(ma a, int k)
{
	if(k<=1)
		return a;
	if(k%2)
		return madd(doit(a,k-1),mpow(a,k));
	ma temp = doit(a,k/2);
	return madd(temp,mmult(temp,mpow(a,k/2)));
}

int main()
{
	int i,j;
	ma a;
	ma s;
	
	e.d[0][0] = e.d[1][1] = 1;
	e.d[0][1] = e.d[1][0] = 0;
	scanf("%d%d%d",&n,&k,&m);
	
	for(i=0; i<n; i++)
	{
		for(j=0; j<n; j++)
		{
			scanf("%d",&a.d[i][j]);
			s.d[i][j] = 0;
		}
	}
	
	s = doit(a,k);
	print(s);
	getchar();getchar();
}