POJ 3710:Matrix Power Series

Description

给出矩阵 \(n*n\) 的 矩阵\(A\) , 求 \(A^1+A^2+A^3...+A^k\)

Solution

首先我们设 \(S_n=\sum_{i=1}^{n}A^i\)
容易得到结论 : \(S_{a+b}=S_{a}*A_{b}+S_{b}\)
于是我们可以把 \(k\) 二进制分解 , 拆成每一个 \(S_{2^i}\) 的形式再按上面的结论合并就行了.
\(S_{2^i}\) 也可以用上述结论倍增求出.
注意这样会多算一个单位矩阵 , 最后减去就行了.

#include<bits/stdc++.h>
using namespace std;
template<class T>void gi(T &x){
	int f;char c;
	for(f=1,c=getchar();c<'0'||c>'9';c=getchar())if(c=='-')f=-1;
	for(x=0;c<='9'&&c>='0';c=getchar())x=x*10+(c&15);x*=f;
}
const int N=35;
int n,k,mod;
struct data{int a[35][35];}A,S,ret;
inline data operator *(const data &p,const data &q){
	data ret;
	for(int i=0;i<n;i++)
		for(int j=0;j<n;j++){
			ret.a[i][j]=0;
			for(int k=0;k<n;k++)
				ret.a[i][j]=(ret.a[i][j]+p.a[i][k]*q.a[k][j])%mod;
		}
	return ret;
}
inline data operator +(const data &p,const data &q){
	data ret;
	for(int i=0;i<n;i++)
		for(int j=0;j<n;j++)
			ret.a[i][j]=(p.a[i][j]+q.a[i][j])%mod;
	return ret;
}
int main(){
  freopen("pp.in","r",stdin);
  freopen("pp.out","w",stdout);
  cin>>n>>k>>mod;
  for(int i=0;i<n;i++)
	  for(int j=0;j<n;j++)scanf("%d",&A.a[i][j]);
  for(int i=0;i<n;i++)
	  for(int j=0;j<n;j++)S.a[i][j]=ret.a[i][j]=(i==j);
  while(k){
	  if(k&1)ret=ret*A+S;
	  S=S*A+S,A=A*A,k>>=1;
  }
  for(int i=0;i<n;i++)
	  for(int j=0;j<n;j++)if(i==j)ret.a[i][j]=(ret.a[i][j]-1+mod)%mod;
  for(int i=0;i<n;i++){
	  for(int j=0;j<n;j++)printf("%d ",ret.a[i][j]);
	  puts("");
  }
  return 0;
}