POJ 3233 Matrix Power Series

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

方法1.分治。

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N=33;
int unit[N][N],init[N][N],f[N][N],n,k,mod;

void plus(int a[N][N],int b[N][N])
{
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
			a[i][j]=(a[i][j]+b[i][j])%mod;
}

void mult(int a[N][N],int b[N][N])
{
	int c[N][N];memset(c,0,sizeof(c));
	for(int i=1;i<=n;i++)
		for(int k=1;k<=n;k++)if(a[i][k])
			for(int j=1;j<=n;j++)
				c[i][j]=(c[i][j]+a[i][k]*b[k][j])%mod;
	memcpy(a,c,sizeof(c));
}

void power_mod(int a[N][N],int b)
{
	int c[N][N];memcpy(c,init,sizeof(c));
	while(b)
	{
		if(b&1)mult(c,a);
		mult(a,a);b=b>>1;
	}
	memcpy(a,c,sizeof(c));
}

void dfs(int a[N][N],int x)
{
	if(x==1)
	{
		memcpy(a,unit,sizeof(unit));
		return;
	}
	int son[N][N],b[N][N];
	memset(son,0,sizeof(son));
	memcpy(b,unit,sizeof(unit));
	dfs(son,x>>1);
	if(x&1)
	{
		power_mod(b,x/2+1);
		plus(a,b);
		plus(a,son);
		mult(son,b);
		plus(a,son);
	}
	else
	{
		power_mod(b,x/2);
		plus(a,son);
		mult(son,b);
		plus(a,son);
	}
}
int main()
{
	memset(f,0,sizeof(f));
	memset(init,0,sizeof(init));
	scanf("%d%d%d",&n,&k,&mod);
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)scanf("%d",&unit[i][j]);
		init[i][i]=1;
	}
	dfs(f,k);
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<n;j++)printf("%d ",f[i][j]);
		printf("%d
",f[i][n]);
	}
	return 0;
}

方法2.我们要求的是等比矩阵，网上大佬提供了一种极妙的构造方法。
$egin{vmatrix}A&amp;1\0&amp;1end{vmatrix}的n次方等于egin{vmatrix}A^n &amp;1+A+A^2+……A^{n-1}\0&amp;1end{vmatrix}$,那么我们可以把矩阵扩大，构造一个$(2n)*（2n)$方阵，求方阵的k+1次方，那么右上角的n*n的方阵-1就是答案。

你可能会问，1怎么用矩阵表示呢？
其实，1在矩阵中指单位矩阵，用I或E表示，上面的-1，+1都是指与单位矩阵的运算。

由于第二行恒为0 1，我们可以不用管它。
也就是求$egin{vmatrix}A&amp;1end{vmatrix}*egin{vmatrix}A&amp;1\0&amp;1end{vmatrix}^k$（但由于这样代码有点复杂，所以笔者就不打了）。

代码：

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N=66;
int E[N][N],f[N][N],n,m,k,mod;
void mult(int a[N][N],int b[N][N])
{
	int c[N][N];memset(c,0,sizeof(c));
	for(int i=1;i<=m;i++)
		for(int k=1;k<=m;k++)if(a[i][k])
			for(int j=1;j<=m;j++)
				c[i][j]=(c[i][j]+a[i][k]*b[k][j])%mod;
	memcpy(a,c,sizeof(c));
}
void power_mod(int a[N][N],int b)
{
	int c[N][N];memcpy(c,E,sizeof(c));
	while(b>0)
	{
		if(b&1)mult(c,a);
		mult(a,a);b=b>>1;
	}
	memcpy(a,c,sizeof(c));
}
int main()
{
	memset(f,0,sizeof(f));
	scanf("%d%d%d",&n,&k,&mod);m=n<<1;
	for(int i=1;i<=m;i++)E[i][i]=1;
	for(int i=1;i<=n;i++)//左上 
		for(int j=1;j<=n;j++)
			scanf("%d",&f[i][j]);
	for(int i=1;i<=n;i++)//右上 
		for(int j=n+1;j<=m;j++)
			f[i][j]=E[i][j-n];
	for(int i=n+1;i<=m;i++)//右下
		for(int j=n+1;j<=m;j++)
			f[i][j]=E[i-n][j-n];
	power_mod(f,k+1);
	for(int i=1;i<=n;i++)
	{
		for(int j=n+1;j<=m;j++)
		{
			f[i][j]-=E[i][j-n];
			if(f[i][j]<0)f[i][j]+=mod;
			if(j==m){printf("%d
",f[i][m]);break;}
			printf("%d ",f[i][j]);
		}
	}
	return 0;
}