【题解】【模板】矩阵级数

【题解】矩阵乘法

给定\(n\)阶矩阵，求\(\Sigma_{i=1}^{k} A^i\)，对\(mod\)取模

此题可以直接分治，但我用纯数学办法做出来了QAQ

在我的另一篇博客里，有这道题的分治做法。

我们直接对题目要求什么设未知数吧,\(sum_i=\Sigma_{j=0}^{i} A^j\),联系无穷级数的套路，从\(0\)开始，可能有比较好的性质(实际上是为了初始化方便)。构造一个矩阵把它们联系到一起。

\[\left[\begin{matrix} A^i\\ S_i\\ \end{matrix}\right] \]

考虑增广

\[\left[\begin{matrix} ? \end{matrix}\right] \times\left[\begin{matrix} A^i\\ S_i\\ \end{matrix}\right]=\left[\begin{matrix} A^{i+1}\\ S_{i+1}\\ \end{matrix}\right] \]

根据矩阵乘法法则，直接待定系数法/观察法解出

\[\left[\begin{matrix} ? \end{matrix}\right]=\left[\begin{matrix} A&0\\ I&I\\ \end{matrix}\right] \]

其中\(\left[\begin{matrix} I \end{matrix}\right]\)是纯量阵。

那么

\[{\left[\begin{matrix} A&0\\ I&I\\ \end{matrix}\right]}\times {\left[\begin{matrix} A^i\\ S_i\\ \end{matrix}\right]} = {\left[\begin{matrix} A^{i+1}\\ S_{i+1}\\ \end{matrix}\right]} \]

所以

\[{\left[\begin{matrix} A&0\\ I&I\\ \end{matrix}\right]}^{k+1}\times {\left[\begin{matrix} I\\ 0\\ \end{matrix}\right]} = {\left[\begin{matrix} A^{k+1}\\ S_{k+1}\\ \end{matrix}\right]} \]

由于

\[S_n=\Sigma_{i=0}^{n}A^i \]

所以

\[\Sigma_{i=1}^{k} A^i=S_{k+1}-I \]

对于

\[{\left[\begin{matrix} A&0\\ I&I\\ \end{matrix}\right]}^{k+1} \]

考虑矩阵快速幂，而最后的答案是

\[{\left[\begin{matrix} S_{k+1} \end{matrix}\right]} \]

直接把

\[{\left[\begin{matrix} A^{k+1}\\ S_{k+1}\\ \end{matrix}\right]} \]

下面的\(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
    
*/