题意:给你一个
n
×
n
n\times n
n×n的矩阵,求它的
k
k
k次方模
1000000007
1000000007
1000000007。
n
≤
50
,
k
≤
2
10000
n\le50,k\le2^{10000}
n≤50,k≤210000。
首先,特征多项式的定义是什么?
我们定义矩阵
A
A
A的特征多项式
f
(
x
)
=
d
e
t
(
A
−
x
I
)
f(x)=det(A-xI)
f(x)=det(A−xI),其中
I
I
I表示单位矩阵,
d
e
t
det
det表示行列式求值。
怎么求出这个特征多项式呢?
我们把
x
=
0..
n
x=0..n
x=0..n全带进去,得到
n
+
1
n+1
n+1个点值,然后再插个值就好了。这一步的时间复杂度是
n
4
n^4
n4。
求出这个特征多项式有什么用?
根据Cayley-Hamilton定理得
f
(
A
)
=
0
f(A)=0
f(A)=0。
因此我们可以让
A
k
A^k
Ak减去若干个
f
(
A
)
f(A)
f(A)。
g
(
x
)
=
x
k
g(x)=x^k
g(x)=xk,
g
(
x
)
g(x)
g(x)可以减去若干个
f
(
x
)
f(x)
f(x),因此
g
(
x
)
=
(
g
m
o
d
f
)
(
x
)
g(x)=(g\ mod\ f)(x)
g(x)=(g mod f)(x)。
所以我们要做的是多项式快速幂+多项式取模。按照快速幂的过程,求出
x
k
x^k
xk的值,在算的过程中对刚刚求出的
f
f
f取模。这一步的时间复杂度为
O
(
n
2
l
o
g
2
k
)
O(n^2log_{2}k)
O(n2log2k)。
总时间复杂度
O
(
n
4
+
n
2
l
o
g
2
k
)
O(n^4+n^2log_{2}k)
O(n4+n2log2k)。
代码
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
const int N=55,mod=1000000007;
int n,len,x[N],y[N];
char k[10005];
int fastpow(int a,int x){
int res=1;
while(x){
if(x&1){
res=1LL*res*a%mod;
}
x>>=1;
a=1LL*a*a%mod;
}
return res;
}
int getinv(int x){
if(x==1){
return 1;
}else{
return fastpow(x,mod-2);
}
}
struct matrix{
int a[N][N];
matrix(){memset(a,0,sizeof(a));}
matrix operator + (const matrix &b) const{
matrix c;
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++){
c.a[i][j]=(a[i][j]+b.a[i][j])%mod;
}
}
return c;
}
matrix operator * (const matrix &b) const{
matrix c;
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++){
for(int k=1;k<=n;k++){
c.a[i][k]=(c.a[i][k]+1LL*a[i][j]*b.a[j][k])%mod;
}
}
}
return c;
}
matrix operator * (const int &x) const{
matrix c=*this;
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++){
c.a[i][j]=1LL*c.a[i][j]*x%mod;
}
}
return c;
}
}a,b;
int det(matrix x){
int res=1;
for(int i=1;i<=n;i++){
if(!x.a[i][i]){
int j=i;
while(j<=n&&!x.a[j][i]){
j++;
}
if(j>n){
return 0;
}
res=-res;
for(int k=1;k<=n;k++){
swap(x.a[i][k],x.a[j][k]);
}
}
int inv=getinv(x.a[i][i]),mul;
for(int j=i+1;j<=n;j++){
mul=1LL*x.a[j][i]*inv%mod;
for(int k=1;k<=n;k++){
x.a[j][k]=(x.a[j][k]-1LL*mul*x.a[i][k]%mod+mod)%mod;
}
}
}
if(res<0){
res+=mod;
}
for(int i=1;i<=n;i++){
res=1LL*res*x.a[i][i]%mod;
}
return res;
}
struct poly{
int len,a[N*2];
poly(){len=0;memset(a,0,sizeof(a));}
void clear(){
while(len&&!a[len]){
len--;
}
}
poly operator + (const poly &b) const{
poly c;
c.len=max(len,b.len);
for(int i=0;i<=c.len;i++){
c.a[i]=(a[i]+b.a[i])%mod;
}
c.clear();
return c;
}
poly operator - (const poly &b) const{
poly c;
c.len=max(len,b.len);
for(int i=0;i<=c.len;i++){
c.a[i]=(a[i]-b.a[i]+mod)%mod;
}
c.clear();
return c;
}
poly operator * (const poly &b) const{
poly c;
c.len=len+b.len;
for(int i=0;i<=len;i++){
for(int j=0;j<=b.len;j++){
c.a[i+j]=(c.a[i+j]+1LL*a[i]*b.a[j])%mod;
}
}
c.clear();
return c;
}
poly operator * (const int &x) const{
poly c=*this;
for(int i=0;i<=c.len;i++){
c.a[i]=1LL*c.a[i]*x%mod;
}
c.clear();
return c;
}
poly operator / (const poly &b) const{
if(len<b.len){
return poly();
}
poly a=*this,c;
c.len=a.len-b.len;
int inv=getinv(b.a[b.len]);
for(int i=c.len;i>=0;i--){
c.a[i]=1LL*a.a[i+b.len]*inv%mod;
for(int j=i+b.len;j>=i;j--){
a.a[j]=(a.a[j]-1LL*c.a[i]*b.a[j-i]%mod+mod)%mod;
}
}
c.clear();
return c;
}
poly operator % (const poly &b) const{
if(len<b.len){
return *this;
}
poly a=*this,c;
c.len=a.len-b.len;
int inv=getinv(b.a[b.len]);
for(int i=c.len;i>=0;i--){
c.a[i]=1LL*a.a[i+b.len]*inv%mod;
for(int j=i+b.len;j>=i;j--){
a.a[j]=(a.a[j]-1LL*c.a[i]*b.a[j-i]%mod+mod)%mod;
}
}
a.clear();
return a;
}
}res;
poly lagrange(){
poly prod,tmp,divi,res;
prod.a[0]=1;
for(int i=0;i<=n;i++){
tmp.len=1;
tmp.a[0]=x[i];
tmp.a[1]=mod-1;
prod=prod*tmp;
}
for(int i=0;i<=n;i++){
tmp.len=1;
tmp.a[0]=x[i];
tmp.a[1]=mod-1;
divi=prod/tmp;
int val=1;
for(int j=0;j<=n;j++){
if(j!=i){
val=1LL*val*(x[j]-x[i]+mod)%mod;
}
}
val=1LL*y[i]*getinv(val)%mod;
res=res+divi*val;
}
return res;
}
int main(){
scanf("%s%d",k+1,&n);
len=strlen(k+1);
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++){
scanf("%d",&a.a[i][j]);
}
}
b=a;
for(int i=0;i<=n;i++){
x[i]=i;
y[i]=det(b);
for(int j=1;j<=n;j++){
b.a[j][j]=(b.a[j][j]-1+mod)%mod;
}
}
res=lagrange();
poly base,prod;
base.len=1;
base.a[1]=1;
prod.a[0]=1;
for(int i=len;i>=1;i--){
if(k[i]>'0'){
prod=prod*base%res;
}
base=base*base%res;
}
matrix mat,ans;
for(int i=1;i<=n;i++){
mat.a[i][i]=1;
}
for(int i=0;i<=prod.len;i++){
ans=ans+mat*prod.a[i];
mat=mat*a;
}
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++){
printf("%d ",ans.a[i][j]);
}
puts("");
}
return 0;
}