POJ3233 - Matrix Power Series - 二分矩阵快速幂

#include<stdio.h>
#include<string.h>
int n,m;
struct JZ
{
	long long jz[50][50];
}s,map;
JZ add(JZ a,JZ b)
{
	JZ ans;
	int i,j;
	memset(ans.jz,0,sizeof(ans.jz));
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			ans.jz[i][j]+=a.jz[i][j]+b.jz[i][j];
		}
		ans.jz[i][j]=ans.jz[i][j]%m;
	}
	return ans;
}
JZ mul(JZ a,JZ b)
{
	JZ ans;
	int i,j,k;
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			ans.jz[i][j]=0;
			for(k=0;k<n;k++)
			{
				ans.jz[i][j]+=a.jz[i][k]*b.jz[k][j];
			}
			ans.jz[i][j]=ans.jz[i][j]%m;
		}
	}
	return ans;
}
JZ pow(long long d)
{
	JZ ans,item;
	item=map;
	int i,j;
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			if(j==i)
			{
				ans.jz[i][j]=1;
			}
			else 
			{
				ans.jz[i][j]=0;
			}
		}
	}
	while(d)
	{
		if(d%2==1)
		{
			ans=mul(ans,item);
		}
		item=mul(item,item);
		d=d/2;
	}
	return ans;
}
JZ half(long long k)
{
	if(k==1)
	{
		return map;
	}
	JZ item,ans;
	item=half(k/2);
	if(k%2==1)//(1+A^(k/2+1))*half(k/2)+A^(k/2+1)
	{
		ans=pow(k/2+1);
		item=add(item,mul(item,ans));
		item=add(item,ans);
	}
	else//(1+A^(k/2))*half(k/2)
	{
		ans=pow(k/2);
		item=add(item,mul(item,ans));
	}
	return item;
}
int main()
{
	int i,j;
	long long k;
	while(scanf("%d%I64d%d",&n,&k,&m)!=EOF)
	{	
		for(i=0;i<n;i++)
		{
			for(j=0;j<n;j++)
			{
				s.jz[i][j]=0;
			}
		}
		for(i=0;i<n;i++)
		{
			for(j=0;j<n;j++)
			{
				scanf("%I64d",&map.jz[i][j]);
			}
		}
		s=half(k);
		for(i=0;i<n;i++)
		{
			for(j=0;j<n;j++)
			{
				if(j!=0)
				{
					printf(" ");
				}
				printf("%I64d",s.jz[i][j]%m);
			}
			puts("");
		}
	}
}