POJ 3233 Matrix Power Series

Description

Given a n × n matrix A and a positive integer k, find the sum S = A + A² + A³ + … + A^k.

Input

The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 10⁹) and m (m < 10⁴). 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

Source

POJ Monthly--2007.06.03, Huang, Jinsong

~~~~~~~~~~~~~~~~~~~~~~~~~~~

倍增+矩阵乘法~

设sum[k]=a+a^1+a^2+...+a^k。

则k为奇数时，sum[k]=sum[k/2]+sum[k/2]*a^(k/2)+a^k；

为偶数时，sum[k]=sum[k/2]+sum[k/2]*a^(k/2)。

倍增即可，非常好写~

快速幂居然没有加if(v&1)……第二次了。

#include<cstdio>
#include<iostream>
using namespace std;

int n,m,mod;

struct node{
	int a[31][31];
}a;

int read()
{
	int x=0,f=1;char ch=getchar();
	while(ch<'0' || ch>'9') {if(ch=='-') f=-1;ch=getchar();}
	while(ch>='0' && ch<='9') {x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
	return x*f;
}

node operator + (node u,node v)
{
	for(int i=1;i<=n;i++)
	  for(int j=1;j<=n;j++) u.a[i][j]=(u.a[i][j]+v.a[i][j])%mod;
	return u;
}

node operator * (node u,node v)
{
	node z;
	for(int i=1;i<=n;i++)
	  for(int j=1;j<=n;j++)
	  {
	  	z.a[i][j]=0;
	  	for(int k=1;k<=n;k++) z.a[i][j]=(z.a[i][j]+u.a[i][k]*v.a[k][j])%mod;
	  }
	return z;
}

node mi(node u,int v)
{
	node z;
	for(int i=1;i<=n;i++)
	  for(int j=1;j<=n;j++) z.a[i][j]=i==j ? 1:0;
	for(;v;v>>=1,u=u*u) if(v&1) z=z*u;
	return z;
}

node solve(int u)
{
	if(u==1) return a;
	node b=solve(u>>1);
	b=b+(b*mi(a,u>>1));
	if(u&1) b=b+mi(a,u);
	return b;
}

int main()
{
	n=read();m=read();mod=read();
	for(int i=1;i<=n;i++)
	  for(int j=1;j<=n;j++) a.a[i][j]=read()%mod;
	a=solve(m);
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<n;j++) printf("%d ",a.a[i][j]);printf("%d\n",a.a[i][n]);
	}
	return 0;
}