• 矩阵倍增 学习总结


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

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

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

    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的

  • 相关阅读:
    重写分词器
    twitterfeed
    What is WSGI?
    lucene 3.3一元切分查询例子
    How to Import Your Blog to Facebook Automaticly Using Notes
    协同过滤 学习笔记
    Python入门练习(一):基于全切分,一元语法模型的汉语分词
    test blog sync to qq microblogging
    solr的用分布式搜索(转)
    盘古分词功能简介
  • 原文地址:https://www.cnblogs.com/joyouth/p/5539655.html
Copyright © 2020-2023  润新知