【题解】矩阵乘法
给定(n)阶矩阵，求(Sigma_{i=1}^{k} A^i)，对(mod)取模
此题可以直接分治，但我用纯数学办法做出来了QAQ
在我的另一篇博客里，有这道题的分治做法。
我们直接对题目要求什么设未知数吧,(sum_i=Sigma_{j=0}^{i} A^j),联系无穷级数的套路，从(0)开始，可能有比较好的性质(实际上是为了初始化方便)。构造一个矩阵把它们联系到一起。
[left[egin{matrix}
A^i\
S_i\
end{matrix}
ight]
]
考虑增广
[left[egin{matrix}
?
end{matrix}
ight] imesleft[egin{matrix}
A^i\
S_i\
end{matrix}
ight]=left[egin{matrix}
A^{i+1}\
S_{i+1}\
end{matrix}
ight]
]
根据矩阵乘法法则，直接待定系数法/观察法解出
[left[egin{matrix}
?
end{matrix}
ight]=left[egin{matrix}
A&0\
I&I\
end{matrix}
ight]
]
其中(left[egin{matrix} I end{matrix} ight])是纯量阵。
那么
[{left[egin{matrix}
A&0\
I&I\
end{matrix}
ight]} imes {left[egin{matrix}
A^i\
S_i\
end{matrix}
ight]} = {left[egin{matrix}
A^{i+1}\
S_{i+1}\
end{matrix}
ight]}
]
所以
[{left[egin{matrix}
A&0\
I&I\
end{matrix}
ight]}^{k+1} imes {left[egin{matrix}
I\
0\
end{matrix}
ight]} = {left[egin{matrix}
A^{k+1}\
S_{k+1}\
end{matrix}
ight]}
]
由于
[S_n=Sigma_{i=0}^{n}A^i
]
所以
[Sigma_{i=1}^{k} A^i=S_{k+1}-I
]
对于
[{left[egin{matrix}
A&0\
I&I\
end{matrix}
ight]}^{k+1}
]
考虑矩阵快速幂，而最后的答案是
[{left[egin{matrix}
S_{k+1}
end{matrix}
ight]}
]
直接把
[{left[egin{matrix}
A^{k+1}\
S_{k+1}\
end{matrix}
ight]}
]
下面的(S_{k+1})输出出来就好了。记得(A)和(S)都是矩阵，所以可以直接输出
考场代码
#include<bits/stdc++.h>
using namespace std;typedef long long ll;
#define DRP(t,a,b) for(register int t=(a),edd=(b);t>=edd;--t)
#define RP(t,a,b) for(register int t=(a),edd=(b);t<=edd;++t)
#define ERP(t,a) for(register int t=head[a];t;t=e[t].nx)
#define midd register int mid=(l+r)>>1
#define TMP template < class ccf >
TMP inline ccf qr(ccf b){
register char c=getchar();register int q=1;register ccf x=0;
while(c<48||c>57)q=c==45?-1:q,c=getchar();
while(c>=48&&c<=57)x=x*10+c-48,c=getchar();
return q==-1?-x:x;}
TMP inline ccf Max(ccf a,ccf b){return a<b?b:a;}
TMP inline ccf Min(ccf a,ccf b){return a<b?a:b;}
TMP inline ccf Max(ccf a,ccf b,ccf c){return Max(a,Max(b,c));}
TMP inline ccf Min(ccf a,ccf b,ccf c){return Min(a,Min(b,c));}
TMP inline ccf READ(ccf* _arr,int _n){RP(t,1,_n)_arr[t]=qr((ccf)1);}
//----------------------template&IO---------------------------
const int maxn=31<<1|1;
ll n,yyb,k;
struct mtx{
ll data[maxn][maxn];
mtx(){memset(data,0,sizeof data);}
inline ll *operator [](int x){return data[x];}
inline mtx operator *=(mtx x){return *this=(*this)*x;}
inline mtx operator ^=(int x){return *this=(*this)^x;}
inline void unis(){
RP(t,1,n)
data[t][t]=1;
}
inline mtx operator * (mtx x){mtx ret;
RP(t,1,n)RP(i,1,n)RP(k,1,n)
(ret[t][i]+=data[t][k]*x[k][i])%=yyb;
return ret;
}
inline mtx operator ^ (int x){
mtx base=*this,ret;ret.unis();
for(register ll t=x;t;t>>=1,base*=base)
if(t&1)ret*=base;
return ret;
}
}zsy,psj,mona;
int main(){
#ifndef ONLINE_JUDGE
freopen("t1.in","r",stdin);
freopen("t1.out","w",stdout);
#endif
n=qr(1ll);
k=qr(1ll)+1ll;
yyb=qr(1ll);
RP(t,1,n){
RP(i,1,n){
zsy[t][i]=qr(1)%yyb;
}
}
RP(t,1,n)
zsy[t+n][t]=zsy[t+n][t+n]=1;
n<<=1;
zsy^=k;
psj.unis();
n>>=1;
RP(t,1,n) zsy[t+n][t]=(zsy[t+n][t]-1+yyb)%yyb;
RP(t,1,n){
RP(i,1,n){
cout<<zsy[t+n][i]<<' ';
}
cout<<endl;
}
return 0;
}
/*
矩阵快速幂板子题
群论板子题
看一下那个<线性代数>的矩阵分割和线性代数的那个方法就会做了
A 0
I=
I I
*/