矩阵倍增 学习总结

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

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

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

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的

时间： 2024-10-06 00:30:42

矩阵倍增 学习总结的相关文章

UVA 11149 - Power of Matrix(矩阵倍增)

UVA 11149 - Power of Matrix 题目链接 题意:给定一个n*n的矩阵A和k,求∑kiAi 思路:利用倍增去搞,∑kiAi=(1+Ak/2)∑k/2iAi,不断二分即可 代码: #include <cstdio> #include <cstring> const int N = 45; int n, k; struct mat { int v[N][N]; mat() {memset(v, 0, sizeof(v));} mat operator * (mat

UVa 11149 矩阵的幂（矩阵倍增法模板题）

https://vjudge.net/problem/UVA-11149 题意: 输入一个n×n矩阵A,计算A+A^2+A^3+...A^k的值. 思路: 矩阵倍增法. 处理方法如下,一直化简下去直到变成A. 代码如下: 1 Matrix solve(Matrix base,int x) 2 { 3 if(x==1)return base; 4 Matrix temp=solve(base,x/2); 5 Matrix sum=add(temp,multi(pow(base,x/2),temp)

UVA - 11149 Power of Matrix（矩阵倍增）

题意:已知N*N的矩阵A,输出矩阵A + A2 + A3 + . . . + Ak,每个元素只输出最后一个数字. 分析: A + A2 + A3 + . . . + An可整理为下式, 从而可以用log2(n)的复杂度算出结果. 注意:输入时把矩阵A的每个元素对10取余,因为若不处理,会导致k为1的时候结果出错. #include<cstdio> #include<cstring> #include<cstdlib> #include<cctype> #in

cojs QAQ的矩阵 题解报告

题目描述非常的清晰 首先我们考虑(A*B)^m的求法,这个部分可以参考BZOJ 杰杰的女性朋友 我们不难发现(A*B)^m=A*(B*A)^(m-1)*B A*B是n*n的矩阵,而B*A是k*k的矩阵,这样就大大缩小了矩阵的大小 因为矩阵乘法满足结合律,我们先对(B*A)做快速幂,之后乘一下就可以了 之后我们考虑如果没有(i-1)^3的这个系数怎么求G(i)的前缀和 因为矩阵乘法满足分配率,我们利用矩阵倍增((B*A)^0+(B*A)^1……+(B*A)^(m-1))之后乘一下就可以了 之后我们

python ocr(光学文字识别)学习笔记 （二）

参考资料:500 lines or less ocr 其中包括神经网络算法的简单介绍,如果看不懂您需要使用谷歌翻译呢 在这一节内容中,我们将对实现这个系统的算法进行分析 设计feedforward ANN(前馈神经网络,也称bp神经网络)时,我们需要考虑以下因素: 1.激活函数的选用 激活函数是结点输出的决策者.我们这个系统将为每个数字输出一个介于0到1的值,值越接近1意味着ann预测的是绘制的数字,越接近0意味着它被预测不是绘制的数字.因此我们将输出接近0或者1的激活函数.我们还需要一个可微分

深度学习笔记(五)：LSTM

深度学习笔记(一):logistic分类 深度学习笔记(二):简单神经网络,后向传播算法及实现 深度学习笔记(三):激活函数和损失函数 深度学习笔记(四):循环神经网络的概念,结构和代码注释 深度学习笔记(五):LSTM 看到一篇讲LSTM非常清晰的文章,原文来自Understanding LSTM Networks , 译文来自理解LSTM网络 Recurrent Neural Networks 人类并不是每时每刻都从一片空白的大脑开始他们的思考.在你阅读这篇文章时候,你都是基于自己已经拥有的

C数组实现矩阵的转置

直接上代码,在代码中有对矩阵的学习,包括初始化学习以及如何使用等. #include <stdio.h> /** * 给出提示,要求输入数组A * ,通过二维数组,进行数组的转置 * 得出数组B,输出结果 * * 该实例主要是为了进行学习二维数组 * @brief main * @return */ int main(void) { /** * 二维数组的初始化: * 1:分行给二维数组赋值 * static int a[3][4] = {{1,2,3,4},{5,6,7,8},{9,10,1

端口相关知识学习笔记

端口相关知识学习笔记 端口相关知识学习笔记 本周主要精力是放在挂接上,所以知识矩阵的学习回归到根本上,所以这周发的学习笔记是关于计算机端口的相关介绍. 有过一些黑客攻击方面知识的读者都会知道,其实那些所谓的黑客并不是像人们想象那样从天而降,而是实实在在从您的计算机"大门"中自由出入.计算机的" 大门"就是我们平常所说的"端口",它包括计算机的物理端口,如计算机的串口.并口.输入/输出设备以及适配器接口等(这些端口都是可见的),但更多的是不可见的软

Codeforces 781D Axel and Marston in Bitland 矩阵 bitset

原文链接https://www.cnblogs.com/zhouzhendong/p/CF781D.html 题目传送门 - CF781D 题意 有一个 n 个点的图,有 m 条有向边,边有两种类型:0 和 1 . 有一个序列,其构造方案为:初始情况为 0 :对于当前串,将当前串的 0 变成 1 ,1 变成 0 ,接到当前串的后面,得到一个新串:不断进行这个操作. 最终得到一个形如 0110100110010110…… 的无限串. 问从节点 1 出发,依次经过上述串对应的 0/1 边,最多能走多