所謂矩陣倍增，就是考試的時候學習的一種新技巧

從字面上就可以理解，利用倍增思想求得我們所需要的矩陣

理論基礎是 矩陣滿足結合律和分配率

?

UVa 11149

裸題，給定矩陣T，求T^1+……+T^n的矩陣

我們都知道利用倍增我們很容易求出T^i (i=2^k)的矩陣，時間復雜度是O(m^3logn)

我們定義一個矩陣為F(k)=T^1+……+T^(2^k)，定義G(k)=T^(2^k)

不難發現F(k+1)=F(k)+F(k)*G(k),G(k+1)=G(k)*G(k)

之后我們用類似快速冪的方式將n分解，當第k位是1的時候我們計算一下F(k)的貢獻并加入答案

時間復雜度O(m^3logn),常數略大

#include<cstdio> #include<cstring> #include<cstdlib> #include<iostream> #include<algorithm> using namespace std;const int mod=10; int n,k; struct Matrix{int a[42][42];void clear(){memset(a,0,sizeof(a));}void print(){for(int i=1;i<=n;++i){for(int j=1;j<=n;++j){printf("%d",a[i][j]);if(j==n)printf("\n");else printf(" ");}}printf("\n");} }A,B,C,I,ans,base;Matrix operator *(const Matrix &A,const Matrix &B){Matrix C;C.clear();for(int i=1;i<=n;++i){for(int j=1;j<=n;++j){for(int k=1;k<=n;++k){C.a[i][j]=C.a[i][j]+A.a[i][k]*B.a[k][j];}C.a[i][j]%=mod;}}return C; } Matrix operator +(const Matrix &A,const Matrix &B){Matrix C;C.clear();for(int i=1;i<=n;++i){for(int j=1;j<=n;++j){C.a[i][j]=A.a[i][j]+B.a[i][j];if(C.a[i][j]>=mod)C.a[i][j]-=mod;}}return C; } void Solve(int p){B=A;C=A;ans.clear();base=I;while(p){if(p&1){ans=ans+B*base;base=base*C;}B=B+B*C;C=C*C;p>>=1;}return; }int main(){while(scanf("%d%d",&n,&k)==2){if(!n)break;for(int i=1;i<=n;++i){for(int j=1;j<=n;++j){scanf("%d",&A.a[i][j]); A.a[i][j]%=mod; }}I.clear();for(int i=1;i<=n;++i)I.a[i][i]=1;Solve(k);ans.print();}return 0; }

5.25考試題T1

我們也是要求A^0+……+A^(n-1)

當時我的做法是利用等比數列求和和矩陣求逆搞一搞搞出來的

實際上那道題目也可以寫矩陣倍增水過去

#include<cstdio> #include<iostream> #include<algorithm> #include<cstdlib> #include<cstring> using namespace std;typedef long long LL; const int mod=1e9+7; int T,n,k; struct Matrix{LL a[2][2];void clear(){memset(a,0,sizeof(a));} }A,B,C,base,ans,I,Ans;Matrix operator *(const Matrix &A,const Matrix &B){Matrix C;C.clear();for(int i=0;i<2;++i){for(int j=0;j<2;++j){for(int k=0;k<2;++k){C.a[i][j]=C.a[i][j]+A.a[i][k]*B.a[k][j]%mod;if(C.a[i][j]>=mod)C.a[i][j]-=mod;}}}return C; } Matrix operator +(const Matrix &A,const Matrix &B){Matrix C;C.clear();for(int i=0;i<2;++i){for(int j=0;j<2;++j){C.a[i][j]=A.a[i][j]+B.a[i][j];if(C.a[i][j]>=mod)C.a[i][j]-=mod;}}return C; } void Solve(int p){B=A;C=A;base=I;ans=I;while(p){if(p&1){ans=ans+B*base;base=base*C;}B=B+B*C;C=C*C;p>>=1;}return; } Matrix pow_mod(Matrix v,int p){Matrix tmp=I;while(p){if(p&1)tmp=tmp*v;v=v*v;p>>=1;}return tmp; }int main(){scanf("%d",&T);A.a[1][0]=A.a[0][1]=A.a[1][1]=1;I.a[0][0]=I.a[1][1]=1;while(T--){scanf("%d%d",&n,&k);if(n==1){printf("1\n");continue;}n--;Solve(n);ans=pow_mod(ans,k);Ans.clear();Ans.a[0][1]=1;Ans=Ans*ans;printf("%lld\n",Ans.a[0][1]);}return 0; }

BZOJ 杰杰的女性朋友

這道題目在BZOJ雙倍經驗QAQ，還有一道是Graph。。

首先我們把out寫成一個n*k的矩陣，然后把in反過來寫成k*n的矩陣

題目中的一坨東西顯然是一堆向量的線性組合，那么用這樣的矩陣描述可以得到

鄰接矩陣G=out*in,那么路徑恰好為d的條數的矩陣就是G^d

也就是(out*in)^d，我們發現out*in得到的是n*n的矩陣，而in*out得到的是k*k的矩陣

n<=1000，k<=20，所以利用結合律，我們可以把式子寫成out*(in*out)^(d-1)*in

然后因為要求<=d的總和，定義T=in*out

根據分配率可以得到out*(T^1+T^2+……+T^(d-1))*in

中間部分利用矩陣倍增即可，最后注意判斷u==v

#include<cstdio> #include<cstring> #include<iostream> #include<algorithm> #include<cstdlib> using namespace std;typedef long long LL; const int mod=1e9+7; int n,m,K,u,v,d; LL out[1010][22],in[22][1010]; LL tmp[1010][22]; struct Matrix{LL a[22][22];void clear(){memset(a,0,sizeof(a));} }A,B,C,I,base,ans;void read(LL &num){num=0;char ch=getchar();while(ch<'!')ch=getchar();while(ch>='0'&&ch<='9')num=num*10+ch-'0',ch=getchar();num%=mod; } Matrix operator *(const Matrix &A,const Matrix &B){Matrix C;C.clear();for(int i=1;i<=K;++i){for(int j=1;j<=K;++j){for(int k=1;k<=K;++k){C.a[i][j]=C.a[i][j]+A.a[i][k]*B.a[k][j]%mod;if(C.a[i][j]>=mod)C.a[i][j]-=mod;}}}return C; } Matrix operator +(const Matrix &A,const Matrix &B){Matrix C;C.clear();for(int i=1;i<=K;++i){for(int j=1;j<=K;++j){C.a[i][j]=A.a[i][j]+B.a[i][j];if(C.a[i][j]>=mod)C.a[i][j]-=mod;}}return C; } void Get_ans(int p){B=A;C=A;ans=I;base=I;while(p){if(p&1){ans=ans+B*base;base=base*C;}B=B+B*C;C=C*C;p>>=1;}return; }int main(){scanf("%d%d",&n,&K);for(int i=1;i<=n;++i){for(int j=1;j<=K;++j)read(out[i][j]);for(int j=1;j<=K;++j)read(in[j][i]);}for(int i=1;i<=K;++i){for(int j=1;j<=K;++j){for(int k=1;k<=n;++k){A.a[i][j]=A.a[i][j]+in[i][k]*out[k][j]%mod;if(A.a[i][j]>=mod)A.a[i][j]-=mod;}}}for(int i=1;i<=K;++i)I.a[i][i]=1;scanf("%d",&m);while(m--){scanf("%d%d%d",&u,&v,&d);if(d==0){printf("%d\n",(u==v?1:0));continue;}d--;Get_ans(d);for(int i=1;i<=n;++i){for(int j=1;j<=K;++j){tmp[i][j]=0;for(int k=1;k<=K;++k){tmp[i][j]=tmp[i][j]+out[i][k]*ans.a[k][j]%mod;if(tmp[i][j]>=mod)tmp[i][j]-=mod;}}}LL S=0;for(int i=1;i<=K;++i){S=S+tmp[u][i]*in[i][v]%mod;if(S>=mod)S-=mod;}if(u==v)S++;if(S>=mod)S-=mod;printf("%lld\n",S);}return 0; }

BZOJ 林中路徑

題目翻譯過來是求sigma(i^2*T^i)

首先如果沒有i^2我們是會做的，我們不妨先考慮i*T^i的情況

定義A(i)=T^i,B(i)=i*T^i

不難發現B(i) = (i-j)*T^i+j*T^i = ((i-j)*T^(i-j)*T^j)+(j*T^j*T^(i-j))

也就是B(i) = B(i-j)*A(j) + B(j)*A(i-j)

這樣我們定義C(i)=i^2*T^i

同理推導一下可以得到

C(i) = C(i-j)*A(j) + C(j)*A(i-j) +2*B(j)*B(i-j)

然后我們求sigma用矩陣倍增搞一搞就可以了

不小心把ik*kj寫成了ij*jk調了一個多小時QAQ

#include<cstdio> #include<cstring> #include<cstdlib> #include<iostream> #include<algorithm> using namespace std;typedef long long LL; int n,m,k,q,u,v; const int mod=1e9+7; struct Matrix{LL a[110][110];LL b[110][110];LL c[110][110];void clear(){memset(a,0,sizeof(a));memset(b,0,sizeof(b));memset(c,0,sizeof(c));}void print(){for(int i=1;i<=n;++i){for(int j=1;j<=n;++j){printf("%lld ",a[i][j]);}printf("\n");}printf("\n");for(int i=1;i<=n;++i){for(int j=1;j<=n;++j){printf("%lld ",b[i][j]);}printf("\n");}printf("\n");for(int i=1;i<=n;++i){for(int j=1;j<=n;++j){printf("%lld ",c[i][j]);}printf("\n");}printf("\n");} }G,A,B,ans,base;Matrix operator *(const Matrix &A,const Matrix &B){Matrix C;C.clear();for(int i=1;i<=n;++i){for(int j=1;j<=n;++j){for(int k=1;k<=n;++k){C.a[i][j]=C.a[i][j]+A.c[i][k]*B.a[k][j]+A.a[i][k]*B.c[k][j]+2LL*A.b[i][k]*B.b[k][j];C.a[i][j]%=mod;C.b[i][j]=C.b[i][j]+A.c[i][k]*B.b[k][j]+A.b[i][k]*B.c[k][j];C.b[i][j]%=mod;C.c[i][j]=C.c[i][j]+A.c[i][k]*B.c[k][j];C.c[i][j]%=mod;}}}return C; } Matrix operator +(const Matrix &A,const Matrix &B){Matrix C;C.clear();for(int i=1;i<=n;++i){for(int j=1;j<=n;++j){C.c[i][j]=A.c[i][j]+B.c[i][j];if(C.c[i][j]>=mod)C.c[i][j]-=mod;C.b[i][j]=A.b[i][j]+B.b[i][j];if(C.b[i][j]>=mod)C.b[i][j]-=mod;C.a[i][j]=A.a[i][j]+B.a[i][j];if(C.a[i][j]>=mod)C.a[i][j]-=mod;}}return C; } void Solve(int p){for(int i=1;i<=n;++i)ans.c[i][i]=1;for(int i=1;i<=n;++i)base.c[i][i]=1;A=G;B=G;while(p){if(p&1){ans=ans+A*base;base=base*B;}A=A+A*B;B=B*B;p>>=1;}return; }int main(){scanf("%d%d%d%d",&n,&m,&k,&q);for(int i=1;i<=m;++i){scanf("%d%d",&u,&v);G.a[u][v]++;G.b[u][v]++;G.c[u][v]++;}Solve(k);while(q--){scanf("%d%d",&u,&v);printf("%d\n",ans.a[u][v]);}return 0; }

BZOJ 4386

由于邊權只有1,2,3，所以可以把每個點弄成3個點，這樣邊權就都是1了，做法跟SCOI 迷路類似

之后我們考慮二分答案，那么我們利用矩陣倍增就可以判斷是否有k個

這樣時間復雜度O((3*n)^3log^2k)顯然會T

首先這道題目我們為了方便統計，可以采用虛擬節點向每個點連邊，然后這個虛擬點上帶個自環

這樣因為有自環的存在我們就不用矩陣倍增，直接快速冪就可以了

不難發現每次二分我們都做了一些重復的工作，我們可以考慮利用倍增

預處理出G^(2^k)，然后從大到小嘗試添加看看是否方案數>=k

這樣時間復雜度少了一個log

注意到矩陣乘法的過程中可能會爆掉long long

但是我們只需要比較跟K的大小就可以了，爆掉long long顯然比K大，所以用-1表示這種情況即可

忘記寫成1LL<<L結果T的窩不知所措

注意到答案上界是3*K,極限情況是兩個點之間相互有長度為3的邊

#include<cstdio> #include<cstring> #include<iostream> #include<cstdlib> #include<algorithm> using namespace std;typedef long long LL; int n,m,u,v,w,N,L; LL K,ans; int idx[42][3]; int deg[122]; struct Matrix{LL a[122][122];void clear(){memset(a,0,sizeof(a));} }G,C,B,A[72]; void read(int &num){num=0;char ch=getchar();while(ch<'!')ch=getchar();while(ch>='0'&&ch<='9')num=num*10+ch-'0',ch=getchar(); } Matrix operator *(const Matrix &A,const Matrix &B){Matrix C;C.clear();for(int i=0;i<=N;++i){for(int j=0;j<=N;++j){for(int k=0;k<=N;++k){if(A.a[i][k]==-1||B.a[k][j]==-1){C.a[i][j]=-1;break;}if(A.a[i][k]==0||B.a[k][j]==0)continue;if(A.a[i][k]>K/B.a[k][j]){C.a[i][j]=-1;break;}C.a[i][j]+=A.a[i][k]*B.a[k][j];if(C.a[i][j]>K){C.a[i][j]=-1;break;} }}}return C; } bool judge(){LL tot=0;for(int i=0;i<=N;++i){if(!deg[i])continue;if(B.a[0][i]==-1)return false;if(B.a[0][i]>K/deg[i])return false;tot+=B.a[0][i]*deg[i];if(tot>=K)return false;}return true; } int main(){read(n);read(m);scanf("%lld",&K);for(int i=1;i<=n;++i){for(int j=0;j<3;++j)idx[i][j]=++N;}for(int i=1;i<=n;++i){for(int j=1;j<3;++j){u=idx[i][j-1];v=idx[i][j];G.a[u][v]++;}G.a[0][idx[i][0]]++;}G.a[0][0]++;for(int i=1;i<=m;++i){read(u);read(v);read(w);--w;G.a[idx[u][w]][idx[v][0]]++;deg[idx[u][w]]++;}for(L=0;(1LL<<L)<=K*3;L++);A[0]=G;ans=1;for(int i=1;i<L;++i)A[i]=A[i-1]*A[i-1];for(int i=0;i<=N;++i)C.a[i][i]=1;for(int i=L-1;i>=0;--i){B=C*A[i];if(judge()){C=B;ans+=(1LL<<i);}}if(ans>K*3)printf("-1\n");else printf("%lld\n",ans);return 0; }

?

倍增是一種思想，利用的是任意自然數都可以被寫成logn位二進制

如果每一位對于答案的貢獻是可以獨立計算的，那么我們倍增預處理每一位的情況之后合并就可以了

通常也可以用于對二分的轉化，求第k大之類的方式

常見的倍增有快速冪，LCA，RMQ，矩陣倍增等等

時間復雜度是log的

轉載于:https://www.cnblogs.com/joyouth/p/5539655.html

與50位技術專家面對面20年技術見證，附贈技術全景圖

總結

以上是生活随笔為你收集整理的矩阵倍增 学习总结的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。