首页 > 代码库 > 矩阵快速幂的一份小结
矩阵快速幂的一份小结
矩阵真是个好东西!虽然矩乘的复杂度有点难看... ...
这几天也做了不少矩阵题目,还是有几道好题目的。不过我打算从入门开始。
- 矩阵乘法:A[i][k]*B[k][j]=C[i][j];(A的第i行的每项依次乘以B的第j列的每项的和)
很显然这是一个n^3的算法,还是比较难看的。
代码就差不多是这样了。
struct Matrix{int T[51][51];}; Matrix Mul(Matrix a,Matrix b,int I,int K,int J) { Matrix S=S0; for(int i=1;i<=I;++i) for(int j=1;j<=J;++j) for(int k=1;k<=K;++k) S.T[i][j]=(S.T[i][j]+a.T[i][k]*b.T[k][j]%Mod)%Mod; return S; }
然后我们发现矩乘有一些性质:
- 不满足交换律。
- 满足结合律。
所以矩乘是可以用来快速幂的!
- 一些题目
题型1.手推矩阵+快速幂
这类题目很多,简单的有求斐波那契数列第多少多少项,难的也有,一般难题重点不在矩阵上。
举几个栗子:
1.斐波那契数列
刚才奶到的求第n项的问题。初始矩阵为
1 | 0 |
0 | 0 |
转移矩阵为
1 | 1 |
1 | 0 |
代码就大概长这样。矩乘的基本操作了。
#include <iostream> #include <cstdio> #include <cstdlib> #include <algorithm> #include <vector> #include <cstring> #include <queue> #define LL long long int #define ls (x << 1) #define rs (x << 1 | 1) #define MID int mid=(l+r)>>1 using namespace std; int n,Mod,ans[2][3][3],S[3][3]; int gi() { int x=0,res=1;char ch=getchar(); while(ch>‘9‘||ch<‘0‘){if(ch==‘-‘)res*=-1;ch=getchar();} while(ch<=‘9‘&&ch>=‘0‘)x=x*10+ch-48,ch=getchar(); return x*res; } void calc(int a,int b) { for(int i=1;i<=2;++i) for(int j=1;j<=2;++j) for(int k=1;k<=2;++k) S[i][j]=(S[i][j]+ans[a][i][k]*ans[b][k][j])%Mod; for(int i=1;i<=2;++i) for(int j=1;j<=2;++j) ans[a][i][j]=S[i][j],S[i][j]=0; } int Qpow(int z) { for(;z;z>>=1) { if(z&1)calc(0,1); calc(1,1); } return ans[0][1][1]; } int main() { int T=gi(); while(T--) { n=gi();Mod=gi(); if(n==0){printf("0\n");continue;} ans[0][1][1]=1;ans[0][1][2]=0;ans[0][2][1]=0;ans[0][2][2]=0; ans[1][1][1]=1;ans[1][1][2]=1;ans[1][2][1]=1;ans[1][2][2]=0; S[1][1]=S[1][2]=S[2][1]=S[2][2]=0; printf("%d\n",Qpow(n)); } return 0; }
2.Xn数列
题目大意就是给你一个A[n+1]=pA[n]+q的简单数列,求(第n项mod m)mod g;
构造等比数列什么的出去,泥萌已经过时了。
手玩方程构造出矩阵。初始矩阵是2*2的
x0 | C |
0 | 0 |
转移矩阵也是2*2的
a | 0 |
1 | 1 |
剩下的也就那样了。我也不知道为什么是大师的题目。
哦对了矩阵乘法要用龟速乘,不然会爆。
#include <iostream> #include <cstdio> #include <cstdlib> #include <algorithm> #include <vector> #include <cstring> #include <queue> #define int long long #define ls (x << 1) #define rs (x << 1 | 1) #define MID int mid=(l+r)>>1 using namespace std; int m,a,c,x0,n,g; int A[5][5],S[5][5],ans[5][5][5]; int gi() { int x=0,res=1;char ch=getchar(); while(ch>‘9‘||ch<‘0‘){if(ch==‘-‘)res*=-1;ch=getchar();} while(ch<=‘9‘&&ch>=‘0‘)x=x*10+ch-48,ch=getchar(); return x*res; } int aXb(int a,int b,int mod) { a%=mod;b%=mod;int ans=0; for(;b;b>>=1) { if(b&1)ans=(ans+a)%mod; a=(a<<1)%mod; } return ans; } void calc(int a,int b) { for(int i=1;i<=2;++i) for(int j=1;j<=2;++j) for(int k=1;k<=2;++k) S[i][j]=(S[i][j]+aXb(ans[a][i][k],ans[b][k][j],m))%m; for(int i=1;i<=2;++i) for(int j=1;j<=2;++j) ans[a][i][j]=S[i][j],S[i][j]=0; } void Qpow(int z) { for(;z;z>>=1) { if(z&1)calc(0,1); calc(1,1); } } void work() { int Ans=ans[0][1][1]%g; cout<<Ans<<endl; } main() { m=gi();a=gi();c=gi();x0=gi();n=gi();g=gi(); ans[0][1][1]=x0;ans[0][1][2]=c; ans[1][1][1]=a;ans[1][2][1]=1;ans[1][2][2]=1; Qpow(n);work();return 0; }
2.要根据输入处理转移矩阵。
1.沼泽鳄鱼
来我们理理思路。
最暴力的暴力都会打吧,随便怎么样的DFS或者BFS。或者你还可以加个记忆化。
然后是一些特殊的点。在没有任何一只鳄鱼的情况下我们就可以用矩阵乘法——初始矩阵是1×n的,(1,ST)=1,其他点都是0的矩阵。转移矩阵就是n*n对应的邻接矩阵。(手玩体会感受一下它的转移)。
下面就是对鳄鱼的处理了。鳄鱼在那个点的话,就不能走到嘛。所以这个转移矩阵的那个对应的"1"就要变成"0"。处理完这一秒之后这样就变成了这一步的转移矩阵。这样每次都在邻接矩阵基础上修改转移矩阵,也能得到正确的答案。但我们没有控制1s的能力,这样会T。
再看看,发现鱼的线路是有规律的,有周期的,而且只有3,4,5三种周期。所以12s后的转移矩阵和12s前的是相同的——一个循环结出现了!于是我们兴高采烈地用上矩阵乘法结合律,把一个循环结乘成一个矩阵,快速幂就好了。至于剩下的一点点余数,暴力就可以啦,反正也就那么点点东西。
#include <iostream> #include <cstdio> #include <cstdlib> #include <algorithm> #define LL long long int using namespace std; struct Matrix{int T[51][51];}S0,Tx[20],Alt,Ans; int n,m,st,ed,KKK,Fish,Map[20][99],map[99][99]; int gi() { int x=0,res=1;char ch=getchar(); while(ch>‘9‘||ch<‘0‘){if(ch==‘-‘)res*=-1;ch=getchar();} while(ch<=‘9‘&&ch>=‘0‘)x=x*10+ch-48,ch=getchar(); return x*res; } Matrix Mul(Matrix a,Matrix b,int I,int K,int J) { Matrix S=S0; for(int i=1;i<=I;++i) for(int j=1;j<=J;++j) for(int k=1;k<=K;++k) S.T[i][j]=(S.T[i][j]+a.T[i][k]*b.T[k][j]%10000)%10000; return S; } Matrix Qpow(int z) { Matrix d=Ans,f=Alt; for(;z;z>>=1,f=Mul(f,f,n,n,n)) if(z&1)d=Mul(d,f,1,n,n); return d; } int main() { n=gi();m=gi();st=gi()+1;ed=gi()+1;KKK=gi(); for(int i=1;i<=m;++i) { int u=gi()+1,v=gi()+1; map[u][v]=map[v][u]=1; } Fish=gi(); for(int i=1;i<=Fish;++i) { int T=gi(); for(int j=0;j<T;++j) { int x=gi()+1; for(int k=j;k<12;k+=T) Map[k][x]=1; } } for(int t=0;t<12;++t) for(int i=1;i<=n;++i) for(int j=1;j<=n;++j) if(map[i][j] && !Map[(t+1)%12][j]) Tx[t].T[i][j]=1; Alt=Tx[0];Ans.T[1][st]=1; for(int i=1;i<12;++i) Alt=Mul(Alt,Tx[i],n,n,n); Ans=Qpow(KKK/12); for(int t=0,s=KKK%12;t<s;++t) Ans=Mul(Ans,Tx[t%12],1,n,n); printf("%d",Ans.T[1][ed]); return 0; }
2.GT考试(光头考试)
题面就是要你找有多少种没出现过一个数字串的数字串。
方案,取膜,匹配,想到DP。
设状态。想知道有没有出现过,就要知道有没有匹配到m位。所以状态就是[i,j],表示主串匹配到第i位,小串匹配到第j位。答案就是sigma F[n,i],i from 1 to m-1;
想转移。F[i,j]表示的既然是最多能匹配到的位置,那就必须用到KMP来求匹配的最长前缀。先预处理出来。转移就明显了,把上一次j的加到这一次相应的j‘里面就好。这么做固然是正确的,但是会T,N<=10^9... ...
找优化。我们发现没一次的转移都可以写成f[i+1,J]=f[i,j1]*...+f[i,j2]*...+...——一阶行列式!矩阵乘法!快速幂!(想想鳄鱼沼泽那题,没有鳄鱼的情况就是一个一样的一阶行列式)。快速幂码完走人。
总结:这题思维和代码都是有的,还要用到KMP这种东西。B站上一次过,在校内OJ被某无良WG卡了输入... ...
#include <iostream> #include <cstdio> #include <cstdlib> #include <algorithm> #include <vector> #include <cstring> #include <queue> #define LL long long int #define ls (x << 1) #define rs (x << 1 | 1) using namespace std; struct Matrix{int T[22][22];}S0,ST,Inv; int N,M,Mod,Next[25],Hash[200],A[200],Ans; int gi() { int x=0,res=1;char ch=getchar(); while(ch>‘9‘||ch<‘0‘){if(ch==‘-‘)res*=-1;ch=getchar();} while(ch<=‘9‘&&ch>=‘0‘)x=x*10+ch-48,ch=getchar(); return x*res; } LL gl() { LL x=0,res=1;char ch=getchar(); while(ch>‘9‘||ch<‘0‘){if(ch==‘-‘)res*=-1;ch=getchar();} while(ch<=‘9‘&&ch>=‘0‘)x=x*10+ch-48,ch=getchar(); return x*res; } Matrix Mul(Matrix a,Matrix b,int I,int K,int J) { Matrix S=S0; for(int i=0;i<I;++i) for(int j=0;j<J;++j) for(int k=0;k<K;++k) S.T[i][j]=(S.T[i][j]+a.T[i][k]*b.T[k][j]%Mod)%Mod; return S; } Matrix Qpow(Matrix d,int z,int len) { Matrix ans=S0; for(int i=0;i<len;++i)ans.T[i][i]=1; for(;z;z>>=1,d=Mul(d,d,len,len,len)) if(z&1)ans=Mul(ans,d,len,len,len); return ans; } int main() { N=gi();M=gi();Mod=gi(); for(int i=1;i<=M;++i)A[i]=getchar()-‘0‘; Next[1]=0;ST.T[0][0]=1; for(int i=2,j=0;i<=M;++i) { while(j&&A[j+1]!=A[i])j=Next[j]; if(A[j+1]==A[i])++j;Next[i]=j; } for(int i=0;i<M;++i) for(int j=0;j<10;++j) { int x=i; while(x && A[x+1]!=j)x=Next[x]; if(A[x+1]==j)Inv.T[i][x+1]=(Inv.T[i][x+1]+1)%Mod; else Inv.T[i][0]=(Inv.T[i][0]+1)%Mod; } Inv=Qpow(Inv,N,M);ST=Mul(ST,Inv,M,M,M); for(int i=0;i<M;++i)Ans=(Ans+ST.T[0][i])%Mod; printf("%d\n",Ans); return 0; }
上面几题都是很水的题目了(by An***大佬,我知道他不会看见),祝大家统统成神犇!
矩阵快速幂的一份小结