所谓矩阵倍增，就是考试的时候学习的一种新技巧
从字面上就可以理解，利用倍增思想求得我们所需要的矩阵
理论基础是 矩阵满足结合律和分配率
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的