题意
给你一个 n × n n\times n n×n的矩阵 A A A,求 A + A 2 + . . . + A k A+A^2+...+A^k A+A2+...+Ak
运算结果对 m m m取模
n ≤ 30 , m < 1 0 4 , k ≤ 1 0 9 n\le30,m\lt10^4,k\le10^9 n≤30,m<104,k≤109
题解
设 E E E为单位矩阵
方法一 分治
令 f ( k ) = A + A 2 + . . . + A k f(k)=A+A^2+...+A^k f(k)=A+A2+...+Ak
那么
f ( k ) = { f ( k 2 ) ( E + A k 2 ) 2 ∣ k f ( k − 1 2 ) ( E + A k − 1 2 ) + A k 2 ∤ k f(k)= \begin{cases} f(\frac k2)(E+A^{\frac k2})& \text{$2\mid k$} \\ f(\frac{k-1}2)(E+A^{\frac{k-1}2})+A^k& \text{$2\nmid k $} \end{cases} f(k)={f(2k)(E+A2k)f(2k−1)(E+A2k−1)+Ak2∣k2∤k
分治计算 f ( k ) f(k) f(k)即可
时间复杂度 O ( n 3 log 2 k ) O(n^3\log^2k) O(n3log2k)
#include<cstdio>
#include<cstring>
#define fp(i,a,b) for(register int i=a,I=b+1;i<I;++i)
#define fd(i,a,b) for(register int i=a,I=b-1;i>I;--i)
#define file(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)
template<class T>inline bool cmax(T&a,const T&b){return a<b?a=b,1:0;}
template<class T>inline bool cmin(T&a,const T&b){return a>b?a=b,1:0;}
using namespace std;
const int N=30+5,inf=~0u>>1;
typedef int arr[N];
int n,Pow,Mod;
inline int plus(int a,int b){return (a+=b)>=Mod?a-Mod:a;}
struct Matrix{
int a[N][N];
Matrix(int d=0){memset(a,0,sizeof a);if(d)fp(i,1,N-1)a[i][i]=1;}
inline int*operator[](int x){return a[x];}
Matrix operator+(Matrix b){
Matrix c;
fp(i,1,n)fp(j,1,n)
c[i][j]=plus(a[i][j],b[i][j]);
return c;
}
inline Matrix operator*(Matrix b)const{
Matrix c;
fp(k,1,n)fp(i,1,n)if(a[i][k])fp(j,1,n)
c[i][j]=plus(c[i][j],a[i][k]*b[k][j]%Mod);
return c;
}
Matrix operator^(int b){
Matrix x(1),A=*this;
for(;b;b>>=1,A=A*A)if(b&1)x=x*A;
return x;
}
}A,S,E(1),Z;
Matrix f(int k){
if(k==1)return A;
return f(k>>1)*(E+(A^(k>>1)))+(k&1?A^k:Z);
}
int main(){
#ifndef ONLINE_JUDGE
file("s");
#endif
scanf("%d%d%d",&n,&Pow,&Mod);
fp(i,1,n)fp(j,1,n)
scanf("%d",&A[i][j]);
S=f(Pow);
fp(i,1,n){
fp(j,1,n)printf("%d ",S[i][j]);
puts("");
}
return 0;
}
方法二 构造
令
S
2
=
[
A
E
0
E
]
S^2= \begin{bmatrix} A&E\\ 0&E \end{bmatrix}
S2=[A0EE]
可以发现
S
2
=
[
A
2
A
+
E
0
E
]
S
3
=
[
A
3
A
2
+
A
+
E
0
E
]
.
.
.
S
k
+
1
=
[
A
k
+
1
A
k
+
A
k
−
1
+
.
.
.
+
A
+
E
0
E
]
S^2= \begin{bmatrix} A^2&A+E\\ 0&E \end{bmatrix} \ \ \ S^3= \begin{bmatrix} A^3&A^2+A+E\\ 0&E \end{bmatrix}\ ...\ \ S^{k+1}= \begin{bmatrix} A^{k+1}&A^k+A^{k-1}+...+A+E\\ 0&E \end{bmatrix}
S2=[A20A+EE] S3=[A30A2+A+EE] ... Sk+1=[Ak+10Ak+Ak−1+...+A+EE]
然后根据分块矩阵理论可以把
S
S
S扩充成一个
2
n
×
2
n
2n\times 2n
2n×2n的新矩阵做矩阵快速幂
记得最后减去单位矩阵
时间复杂度 O ( n 3 log k ) O(n^3\log k) O(n3logk)
#include<cstdio>
#include<cstring>
#define fp(i,a,b) for(register int i=a,I=b+1;i<I;++i)
#define fd(i,a,b) for(register int i=a,I=b-1;i>I;--i)
#define file(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)
template<class T>inline bool cmax(T&a,const T&b){return a<b?a=b,1:0;}
template<class T>inline bool cmin(T&a,const T&b){return a>b?a=b,1:0;}
using namespace std;
const int N=30*2+5,inf=~0u>>1;
typedef int arr[N];
int n,Pow,Mod;
struct Matrix{
int a[N][N];
Matrix(int d=0){memset(a,0,sizeof a);if(d)fp(i,1,n)a[i][i]=1;}
inline int*operator[](int x){return a[x];}
inline Matrix operator*(Matrix b)const{
Matrix c;
fp(k,1,n)fp(i,1,n)if(a[i][k])fp(j,1,n)
c[i][j]=(c[i][j]+a[i][k]*b[k][j])%Mod;
return c;
}
Matrix operator^(int b){
Matrix x(1),A=*this;
for(;b;b>>=1,A=A*A)if(b&1)x=x*A;
return x;
}
}S;
int main(){
#ifndef ONLINE_JUDGE
file("s");
#endif
scanf("%d%d%d",&n,&Pow,&Mod);
fp(i,1,n)fp(j,1,n){
scanf("%d",&S[i][j]);
if(i==j)S[i+n][j+n]=S[i][j+n]=1;
}
n<<=1;
S=S^(Pow+1);
n>>=1;
fp(i,1,n){
fp(j,1,n)printf("%d ",(S[i][j+n]-(i==j)+Mod)%Mod);
puts("");
}
return 0;
}