几年前刚学矩阵快速幂的时候点开了这道题，写了半天始终TLE，一看题解发现完全看不懂……时隔多年终于A掉了这道题.
因为代码习惯，下文中用$n$表示题目中的$k$，用$s$表示题目中的$n$.
设$M$的特征多项式为$f(\lambda)$，设$\lambda^s=f(\lambda)Q(\lambda)+R(\lambda)$，其中$\deg R<\deg f$，那么因为$f(M)=0$，所以$M^s=R(M)$.
求特征多项式的话，可以把$\lambda=0,1,\ldots,n$带入$\det(M-\lambda I)$，然后插值，计算$R(\lambda)$就直接在模${f(\lambda)}$意义下快速幂就行了.
发现多项式取模的一个比较简单的$\mathcal O(n^2)$写法，设$f(\lambda)=\sum _{i=0}^na _i\lambda^i$，那么$M^n=-\sum _{i=0}^{n-1}\frac{a _i}{a _n}M^i$，利用这个把次数不小于$\deg f$的项全部变成次数小于$\deg f$的项就可以了.
代码写得比较凌乱.
代码：1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
using namespace std;
struct matrix{
ele n,m,a[maxn][maxn];
inline void I(){
memset(a,0,sizeof(a));
for (int i=0; i<n; ++i) a[i][i]=1;
}
};
inline void operator+=(matrix&a,const matrix&b){
for (int i=0; i<a.n; ++i)
for (int j=0; j<a.m; ++j) (a.a[i][j]+=b.a[i][j])%=MOD;
}
inline void operator*=(matrix&a,ele b){
for (int i=0; i<a.n; ++i)
for (int j=0; j<a.m; ++j) a.a[i][j]=(ll)a.a[i][j]*b%MOD;
}
inline void operator*=(matrix&a,const matrix&b){
static matrix c;
memset(c.a,0,sizeof(c.a));
c.n=a.n; c.m=b.m;
for (int i=0; i<a.n; ++i)
for (int j=0; j<a.m; ++j)
for (int k=0; k<b.m; ++k)
(c.a[i][k]+=(ll)a.a[i][j]*b.a[j][k]%MOD)%=MOD;
a=c;
}
char s[maxm];
ele n,f[maxn][maxn],a[maxn][maxn],b[maxn][maxn],S[maxn][maxn],fac[maxn],ifac[maxn],c[maxn],d[maxn],e[maxn];
matrix A,B,C;
inline ele pw(ele a,ele x){
ele ans=1,tmp=a%MOD;
for (; x; x>>=1,tmp=(ll)tmp*tmp%MOD)
if (x&1) ans=(ll)ans*tmp%MOD;
return ans;
}
inline ele det(ele n,ele a[maxn][maxn]){
ele ans=1;
for (int i=0; i<n; ++i)
for (int j=i+1; j<n; ++j)
while (a[j][i]){
ans=MOD-ans;
ele tmp=a[i][i]/a[j][i];
for (int k=i; k<n; ++k){
(a[i][k]+=MOD-(ll)a[j][k]*tmp%MOD)%=MOD;
swap(a[i][k],a[j][k]);
}
}
for (int i=0; i<n; ++i) ans=(ll)ans*a[i][i]%MOD;
return ans;
}
inline void mul(ele *a,ele *b){
static ele t[maxn];
memset(t,0,sizeof(t));
for (int i=0; i<=n; ++i)
for (int j=0; j<=n; ++j) (t[i+j]+=(ll)a[i]*b[j]%MOD)%=MOD;
ele invcn=pw(c[n],MOD-2);
for (int i=n*2; i>=n; --i){
ele tmp=(ll)invcn*t[i]%MOD;
for (int j=1; j<=n; ++j)
(t[i-j]+=MOD-(ll)tmp*c[n-j]%MOD)%=MOD;
t[i]=0;
}
memcpy(a,t,sizeof(ele)*(n+1));
}
int main(){
scanf("%s%d",s,&n);
for (int i=0; i<n; ++i)
for (int j=0; j<n; ++j) scanf("%d",&a[i][j]);
for (int i=0; i<=n; ++i){
memcpy(b,a,sizeof(a));
f[0][i]=det(n,b);
for (int j=0; j<n; ++j) (a[j][j]+=MOD-1)%=MOD;
}
for (int i=1; i<=n; ++i)
for (int j=0; j<=n-i; ++j) f[i][j]=(f[i-1][j+1]-f[i-1][j]+MOD)%MOD;
for (int i=0; i<n; ++i) (a[i][i]+=n+1)%=MOD;
memset(S,0,sizeof(S));
S[0][0]=1;
for (int i=1; i<=n; ++i){
S[i][0]=0;
for (int j=1; j<=i; ++j)
S[i][j]=(S[i-1][j-1]+(ll)S[i-1][j]*(i-1)%MOD)%MOD;
}
fac[0]=1;
for (int i=1; i<=n; ++i) fac[i]=(ll)fac[i-1]*i%MOD;
ifac[n]=pw(fac[n],MOD-2);
for (int i=n-1; ~i; --i) ifac[i]=(ll)ifac[i+1]*(i+1)%MOD;
memset(c,0,sizeof(c));
for (int i=0; i<=n; ++i){
for (int j=i,ct=1; ~j; --j,ct=MOD-ct)
(c[j]+=(ll)S[i][j]*f[i][0]%MOD*ifac[i]%MOD*ct%MOD)%=MOD;
}
memset(d,0,sizeof(d)); d[1]=1;
memset(e,0,sizeof(e)); e[0]=1;
ele L=strlen(s);
for (int i=L-1; ~i; --i){
if (s[i]=='1') mul(e,d);
mul(d,d);
}
A.n=A.m=n;
for (int i=0; i<n; ++i)
for (int j=0; j<n; ++j) A.a[i][j]=a[i][j];
B.n=B.m=n; B.I();
C.n=C.m=n; memset(C.a,0,sizeof(C.a));
for (int i=0; i<n; ++i){
if (e[i]){
B*=e[i];
C+=B;
B*=pw(e[i],MOD-2);
}
B*=A;
}
for (int i=0; i<n; ++i){
for (int j=0; j<n; ++j)
printf(j?" %d":"%d",C.a[i][j]);
puts("");
}
return 0;
}