當前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

矩阵倍增学习总结

發布時間：2024/8/26 编程问答 37 豆豆

生活随笔收集整理的這篇文章主要介紹了矩阵倍增学习总结小編覺得挺不錯的,現在分享給大家,幫大家做個參考.

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

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

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

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年技術見證，附贈技術全景圖