矩阵倍增 学习总结

Posted 2020-07-16 _Vertical

所谓矩阵倍增，就是考试的时候学习的一种新技巧

从字面上就可以理解，利用倍增思想求得我们所需要的矩阵

理论基础是 矩阵满足结合律和分配率

 

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的