所谓矩阵倍增,就是考试的时候学习的一种新技巧
从字面上就可以理解,利用倍增思想求得我们所需要的矩阵
理论基础是 矩阵满足结合律和分配率
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(" "); else printf(" "); } }printf(" "); } }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 ");continue;} n--;Solve(n); ans=pow_mod(ans,k); Ans.clear();Ans.a[0][1]=1; Ans=Ans*ans; printf("%lld ",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 ",(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 ",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(" "); }printf(" "); for(int i=1;i<=n;++i){ for(int j=1;j<=n;++j){ printf("%lld ",b[i][j]); }printf(" "); }printf(" "); for(int i=1;i<=n;++i){ for(int j=1;j<=n;++j){ printf("%lld ",c[i][j]); }printf(" "); }printf(" "); } }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 ",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 "); else printf("%lld ",ans); return 0; }
倍增是一种思想,利用的是任意自然数都可以被写成logn位二进制
如果每一位对于答案的贡献是可以独立计算的,那么我们倍增预处理每一位的情况之后合并就可以了
通常也可以用于对二分的转化,求第k大之类的方式
常见的倍增有快速幂,LCA,RMQ,矩阵倍增等等
时间复杂度是log的