传送门
题意
给出n,m,k,求
[sum_{i=1}^kA^i
]
A是矩阵
分析
我们首先会想到等比公式，然后得到这样一个式子：
[frac{A^{k+1}-E}{A-E}
]
发现要用矩阵除法，可以用求矩阵逆来做，现在我们换一种做法，我们发现有这样一个性质：
[left[
egin{matrix}
A & E \
0 & E \
end{matrix}
ight]^n=
left[
egin{matrix}
A^{n} & sum_{i=0}^{n-1}A^i \
0 & E \
end{matrix}
ight]
]
那么我们将原先矩阵扩大成四倍，对矩阵求k+1次幂，然后取右上角减去一个单位阵即可
trick
代码
#include<cstdio>
#include<cstring>
#include<iostream>
using namespace std;
struct Matrix
{
int matrix[80][80];
}ans,ret;
int t,n,k,mod;
Matrix multi(Matrix a,Matrix b)
{
Matrix tmp;
for(int i=1;i<=n;++i)for(int j=1;j<=n;++j)
{
tmp.matrix[i][j]=0;
for(int k=1;k<=n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
}
for(int i=1;i<=n;++i)for(int j=1;j<=n;++j)
{
for(int k=n+1;k<=2*n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
}
for(int i=1;i<=n;++i)for(int j=n+1;j<=2*n;++j)
{
tmp.matrix[i][j]=0;
for(int k=1;k<=n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
}
for(int i=1;i<=n;++i)for(int j=n+1;j<=2*n;++j)
{
for(int k=n+1;k<=2*n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
}
for(int i=1+n;i<=2*n;++i)for(int j=1;j<=n;++j)
{
tmp.matrix[i][j]=0;
for(int k=1;k<=n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
}
for(int i=1+n;i<=2*n;++i)for(int j=1;j<=n;++j)
{
for(int k=n+1;k<=2*n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
}
for(int i=n+1;i<=2*n;++i)for(int j=n+1;j<=2*n;++j)
{
tmp.matrix[i][j]=0;
for(int k=1;k<=n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
}
for(int i=1+n;i<=2*n;++i)for(int j=n+1;j<=2*n;++j)
{
for(int k=n+1;k<=2*n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
}
return tmp;
}
void fast_mod(int p)
{
memset(ans.matrix,0,sizeof(ans.matrix));
for(int i=1;i<=2*n;++i) ans.matrix[i][i]=1;
for(;p;p>>=1,ret=multi(ret,ret)) if(p&1) ans=multi(ans,ret);
}
int main()
{
//freopen("data.in","w",stdout);
while(scanf("%d %d %d",&n,&k,&mod)!=EOF)
{
for(int i=1;i<=n;++i)for(int j=1;j<=n;++j) scanf("%d",&ret.matrix[i][j]);
for(int i=1;i<=n;++i)for(int j=n+1;j<=2*n;++j) if((j-i)==n) ret.matrix[i][j]=1;else ret.matrix[i][j]=0;
for(int i=n+1;i<=2*n;++i)for(int j=1+n;j<=2*n;++j) if(i==j) ret.matrix[i][j]=1;else ret.matrix[i][j]=0;
for(int i=n+1;i<=2*n;++i)for(int j=1;j<=n;++j) ret.matrix[i][j]=0;
//for(int i=1;i<=2*n;++i)for(int j=1;j<=2*n;++j) printf("%d%c",ret.matrix[i][j],j==2*n?'
':' ');
fast_mod(k+1);
//for(int i=1;i<=2*n;++i)for(int j=1;j<=2*n;++j) printf("%d%c",ans.matrix[i][j],j==2*n?'
':' ');
for(int i=1;i<=n;++i)for(int j=n+1;j<=2*n;++j) if((j-i)==n) (ans.matrix[i][j]+=(mod-1))%=mod;
for(int i=1;i<=n;++i)for(int j=1+n;j<=2*n;++j) printf("%d%c",ans.matrix[i][j],j==2*n?'
':' ');
}
return 0;
}