高级算法指北——矩阵乘法、矩阵快速幂及其应用
I 走近矩阵乘法
i 定义
当矩阵 \(A\) 和矩阵 \(B\) 满足矩阵 \(A\) 的列数等于矩阵 \(B\) 的行数时，两个矩阵可以相乘。设 \(A\) 是一个 \(n\times m\) 的矩阵，\(B\) 是一个 \(m \times k\) 的矩阵，若 \(C=A\times B\)，则 \(C\) 是一个 \(n\times k\) 的矩阵。
具体地，矩阵 \(C\) 满足：
即矩阵 \(C\) 的第 \(i\) 行 \(j\) 列为 \(A\) 的第 \(i\) 行与 \(B\) 的第 \(j\) 列对应位置一一相乘所得值的和。简记为“前行乘后列，对应相乘再求和”。
代码很简单，直接 \(n^3\) 暴力求和即可。我们可以把乘法操作封装在一个结构体里面。
struct matrix{
int a[N][N];
friend matrix operator*(matrix x,matrix y){
matrix res;
rep(i,1,n){
rep(j,1,n)
res.a[i][j]=0;
}
rep(i,1,n){
rep(j,1,n){
rep(k,1,n)
res.a[i][j]+=x.a[i][k]*y.a[k][j],res.a[i][j]%=mo;
}
}
return res;
}
}
特别地，有一类 \(n\times n\) 的矩阵，满足它乘上任何矩阵（或任何矩阵乘上它）后的得到的矩阵不变。我们称之为单位矩阵。单位矩阵除位于左上至右下的对角线上的值为 \(1\) 以外，其余的值均为 \(0\)。
ii 矩阵乘法的运算律与矩阵快速幂
矩阵乘法不满足交换律，由定义中行和列的差异可以明显地看出。
矩阵乘法满足结合律。因此，我们可以使用快速幂来优化矩阵的乘方计算，其代码逻辑与普通的快速幂相同，只需要将初始值设为单位矩阵即可。
matrix quick_power(matrix base,int x){
matrix ret;
rep(i,1,n){//初始化为单位矩阵
rep(j,1,n){
if(i==j)ret.a[i][j]=1;
else ret.a[i][j]=0;
}
}
while(x){
if(x&1)ret=ret*base;
base=base*base;
x>>=1;
}
return ret;
}
II 矩阵快速幂的应用
i 矩阵快速幂优化 dp
对于一类转移来源单一且相似、某一维状态范围极大的 dp，我们可以根据转移式构造方阵，使用矩阵乘法表示其转移过程，并使用矩阵快速幂加速转移即可。
- 例 \(1\): 矩阵加速（数列）
已知一个数列 \(a\)，它满足：
求 \(a\) 数列的第 \(n\) 项对 \(10^9+7\) 取余的值。
【解题】
考虑到每一次转移涉及到前三个 \(a\) 的值，考虑将它们放入一个 \(1\times 3\) 的矩阵中去。
转移时，我们显然要将其变为下面这个矩阵：
即
考虑构造一个 \(3\times 3\) 的方阵，使得我们可以通过让原来的矩阵乘上这个方阵，得到新的矩阵。使用让一个原来的矩阵中的元素与方阵中一列对应，通过待定系数的方式，我们可以轻易求出转移方阵 \(G\) 如下：
此外，不难构造出初始矩阵 \(S_1\)：
则有 \(S_i=S_1\times G^{i-1}\)。用矩阵快速幂即可快速求出 \(S_i\) 的值。对本题而言，求出 \(S_{n-2}\) 即可。
- 例 \(2\): [HNOI2010] 公交线路
小Z所在的城市有N个公交车站，排列在一条长(N-1)km的直线上，从左到右依次编号为1到N，相邻公交车站间的距离均为1km。 作为公交车线路的规划者，小Z调查了市民的需求，决定按下述规则设计线路：
- 设共K辆公交车，则1到K号站作为始发站，N-K+1到N号台作为终点站。
- 每个车站必须被一辆且仅一辆公交车经过（始发站和终点站也算被经过）。
- 公交车只能从编号较小的站台驶往编号较大的站台。
- 一辆公交车经过的相邻两个站台间距离不得超过Pkm。
在最终设计线路之前，小Z想知道有多少种满足要求的方案。由于答案可能很大，你只需求出答案对30031取模的结果。
N<=10^9，1<P<=10，K<N，1<K<=P，K<=8
【解题】
不难想到设 \(dp_{i,j}\) 表示考虑到第 \(i\) 个车站，包含 \(i\) 的后 \(p\) 个公交车站停靠的那个车目前是否在那里停最后一个站的状态为 \(j\) 时的方案总数。初始位置为 \(dp_{k,2^k-1}\)，终止位置为 \(dp_{n,2^k-1}\)。考虑转移时，有 \(dp_{i,nw}=\sum dp_{i-1,st}\)，其中 \(nw\) 状态能由 \(st\) 状态左移 \(1\) 位后移动一个车到最低位而得到。不难发现能转移到 \(nw\) 状态种类是固定的，与 \(i\) 无关，故考虑构造方阵帮助转移。
我们可以将所有状态下的 \(dp\) 值排成一横排，形成一个矩阵。对于方阵的第 \(i\) 行第 \(j\) 列上的数，由于它会和原矩阵中第 \(i\) 个数（原矩阵和新矩阵均只有 \(1\) 行）相乘，贡献累加到新矩阵中第 \(j\) 个位置上去，所以若状态 \(i\) 能转移到状态 \(j\)，则值为 \(1\)，否则为 \(0\)。初始矩阵中仅有 \(2^k-1\) 那个状态对应位置上的值为 \(1\)，最终答案即为乘出来的矩阵上 \(2^k-1\) 那个状态对应位置上的值。
int n,k,p,m;//状态总数为m
struct matrix{
int a[N][N];
friend matrix operator *(matrix x,matrix y){
matrix ret;
rep(i,1,m){
rep(j,1,m)
ret.a[i][j]=0;
}
rep(i,1,m){
rep(j,1,m){
rep(l,1,m){
ret.a[i][j]+=x.a[i][l]*y.a[l][j];
ret.a[i][j]%=mo;
}
}
}
return ret;
}
}mat,pre;
int sit[N],id[(1<<11)+5];//有效状态 为了左移一位后不爆炸，多开一倍
matrix quick_power(matrix base,int x){
matrix ret;
rep(i,1,m){
rep(j,1,m){
if(i!=j)ret.a[i][j]=0;
else ret.a[i][j]=1;
}
}
while(x){
if(x&1)ret=ret*base;
base=base*base;
x>>=1;
}
return ret;
}
signed main(){
read(n),read(k),read(p);
rep(i,1,(1<<p)-1){//统计有效状态
if(!(i&1))continue;
if(__builtin_popcount(i)!=k)continue;
sit[++m]=i,id[i]=m;
}
rep(i,1,m){
rep(j,1,m)
pre.a[i][j]=mat.a[i][j]=0;
}
pre.a[1][1]=1;
rep(i,1,m){//初始时所有1排在最后，是最小的（最靠前的）那个可行状态.
int st=sit[i];
rep(j,1,p){
if(!((st>>(j-1))&1))continue;
int nw=((st^(1<<(j-1)))<<1)|1;
if(!id[nw])continue;
mat.a[i][id[nw]]=1;
}
}
pre=pre*quick_power(mat,n-k);
printf("%lld\n",pre.a[1][1]);//最终状态的所有1排在最后，仍然是最小的那个可行状态.
return 0;
}
ii 邻接矩阵乘法解决图论问题
- 定长路径数统计问题：[ZJOI2005] 沼泽鳄鱼
豆豆先生酷爱冒险，他一听说这个消息，立马赶到了池塘，想做第一个在桥上旅游的人。池塘里有食人鱼，食人鱼的行进路线有周期性，这个周期只可能是 \(2\)、\(3\) 或者 \(4\) 个单位时间。每个单位时间里，食人鱼可以从一个石墩游到另一个石墩。每到一个石墩，如果上面有人它就会实施攻击，否则继续它的周期运动。如果没有到石墩，它是不会攻击人的。
借助先进的仪器，豆豆很快就摸清了所有食人鱼的运动规律，他要开始设计自己的行动路线了。每个单位时间里，他只可以沿着石桥从一个石墩走到另一个石墩，而不可以停在某座石墩上不动，因为站着不动还会有其它危险。如果豆豆和某条食人鱼在同一时刻到达了某座石墩，就会遭到食人鱼的袭击，他当然不希望发生这样的事情。
现在豆豆已经选好了两座石墩 \(\mathrm{Start}\) 和 \(\mathrm{End}\)，他想从 \(\mathrm{Start}\) 出发，经过 \(K\) 个单位时间后恰好站在石墩 \(\mathrm{End}\) 上。假设石墩可以重复经过（包括 \(\mathrm{Start}\) 和 \(\mathrm{End}\)），他想请你帮忙算算，这样的路线共有多少种（当然不能遭到食人鱼的攻击）。
对于 \(100 \%\) 的数据，\(1 \leq N \leq 50\)，\(1 \leq K \leq 2 \times 10^9\)，\(1 \leq \mathrm{NFish} \leq 20\)。
【解题】
首先考虑没有食人鱼的时候怎么做。
对于一个邻接矩阵 \(E\)，\(E_{i,j}=1\) 表示两者之间有边可以直达。同样的，它也可以表示从 \(i\) 出发到 \(j\)，走 \(1\) 条边的方案数。考虑我们怎么求从 \(i\) 出发到 \(j\)，走 \(2\) 条边的方案数，显然可以枚举一个点 \(k\)，方案数为 \(\sum E_{i,k}\times E_{k,j}\)，不难发现这和矩阵乘法的式子一致。因而我们可以直接将邻接矩阵平方，得到代表从 \(i\) 出发到 \(j\)，走 \(2\) 条边的方案数的矩阵 \(E_2\)。进一步地，设表示任意两点间走 \(x\) 条边到达的方案数的矩阵为 \(E_x\)，我们不难发现 \(E_i\times E_j=E_{i+j}\)，也就可以得到 \(E_i=E_1^i\)。这样，我们就可以利用矩阵快速幂来快速求解了。
在有食人鱼的情况下，我们发现食人鱼的游荡周期为 \(12\)。因此，我们可以先对于前 \(12\) 个单位时间中的每一个，处理出食人鱼在的所有的石墩。则那一个时刻的矩阵应当为在原来 \(E_1\) 的基础上，将所有有食人鱼在的石墩对应的列设为 \(0\) 之后得到的矩阵（因为有食人鱼在的石墩，即使有边连接，也到不了）。根据前面说的原理，我们将这 \(12\) 个时间的对应矩阵按顺序依次（注意矩阵乘法不符合交换律！）乘起来，便得到在前 \(12\) 分钟里，从任意一点出发到任意一点，走 \(12\) 条边的方案数矩阵。接下来的时间就是对周期的重复，直接对求出来的一个周期里的方案数矩阵做快速幂即可。最后的不满一个周期的时间可以再暴力乘。
时间复杂度 \(O(N^3\log k)\)。
#include<bits/stdc++.h>
#define rep(i,j,k) for(int i=j;i<=k;i++)
#define repp(i,j,k) for(int i=j;i>=k;i--)
#define ls(x) x*2
#define rs(x) x*2+1
#define mp make_pair
#define fir first
#define sec second
#define pvv pair<vector<int>,vector<int>>
#define lowbit(x) x&-x
using namespace std;
typedef long long ll;
const int N=52,M=1e5+5,S=(1<<17)+2,mo=10000,inf=1e9+7;
void read(int &p){
int x=0,w=1;
char ch=0;
while(!isdigit(ch)){
if(ch=='-')w=-1;
ch=getchar();
}
while(isdigit(ch)){
x=(x<<1)+(x<<3)+ch-'0';
ch=getchar();
}
p=x*w;
}
int n,m,q,st,ed,k;
struct matrix{
int a[N][N];
friend matrix operator*(matrix x,matrix y){
matrix res;
rep(i,1,n){
rep(j,1,n)
res.a[i][j]=0;
}
rep(i,1,n){
rep(j,1,n){
rep(k,1,n)
res.a[i][j]+=x.a[i][k]*y.a[k][j],res.a[i][j]%=mo;
}
}
return res;
}
}e,sit[15],sum;
vector<int>shark[22];
int t[22];
matrix quick_power(matrix base,int x){
matrix ret;
rep(i,1,n){
rep(j,1,n){
if(i==j)ret.a[i][j]=1;
else ret.a[i][j]=0;
}
}
while(x){
if(x&1)ret=ret*base;
base=base*base;
x>>=1;
}
return ret;
}
int main(){
read(n),read(m),read(st),read(ed),read(k),st++,ed++;
rep(i,1,m){
int x,y;
read(x),read(y),x++,y++;
e.a[x][y]=e.a[y][x]=1;
}
read(q);
rep(i,1,q){
read(t[i]);
rep(j,1,t[i]){
int x;
read(x),x++;
shark[i].push_back(x);
}
}
rep(i,1,12){//处理前12分钟的
sit[i]=e;
rep(j,1,q){
int tms=i%t[j],p=shark[j][tms];
rep(l,1,n)
sit[i].a[l][p]=0;
}
}
sum=sit[1];
rep(i,2,12)
sum=sum*sit[i];
matrix ans=quick_power(sum,k/12);
k%=12;
rep(i,1,k)
ans=ans*sit[i];
printf("%d",ans.a[st][ed]);
return 0;
}
- 定边数的最短路问题：[USACO07NOV] Cow Relays G
给定一张 \(T\) 条边的无向连通图，求从 \(S\) 到 \(E\) 经过 \(N\) 条边的最短路长度。
对于所有的数据，保证 \(1\le N\le 10^6\)，\(2\le T\le 100\)。
所有的边保证 \(1\le u,v\le 1000\)，\(1\le w\le 1000\)。
【解题】
与上面的方法类似，由于 floyd 求最短路算法的过程十分类似于矩阵乘法，只需要将乘积加和变为加和求 min。加和求 min 的操作显然也有结合律，因而我们仍然可以使用邻接矩阵做乘法，快速幂优化的方式进行求解，只需要将矩阵乘法转为加和求最小值，并在矩阵中存最短路长度即可。
注意在矩阵快速幂时，初始矩阵应当设成作为底数的那个矩阵，且幂次值同步 \(-1\)；而不能用全为 INF 的矩阵替代。
int st,ed,n,m,k,p[N];
struct matrix{
int a[N][N];
friend matrix operator*(matrix x,matrix y){
matrix res;
rep(i,1,n){
rep(j,1,n)
res.a[i][j]=inf;
}
rep(i,1,n){
rep(j,1,n){
rep(k,1,n)
res.a[i][j]=min(res.a[i][j],x.a[i][k]+y.a[k][j]);
}
}
return res;
}
}dis;
struct edge{
int frm,to,v;
}e[M];
matrix quick_power(matrix base,int x){
matrix ret=base;//注意这里的初始值设定，显然不能为inf.
x--;
while(x){
if(x&1)ret=ret*base;
base=base*base;
x>>=1;
}
return ret;
}
int main(){
read(k),read(m),read(st),read(ed);
rep(i,1,m){
read(e[i].v),read(e[i].frm),read(e[i].to);
p[++n]=e[i].frm,p[++n]=e[i].to;
}
sort(p+1,p+n+1),n=unique(p+1,p+n+1)-p-1;
rep(i,1,n){
rep(j,1,n)
dis.a[i][j]=inf;
}
rep(i,1,m){
int x,y;
x=lower_bound(p+1,p+n+1,e[i].frm)-p,y=lower_bound(p+1,p+n+1,e[i].to)-p;
dis.a[x][y]=dis.a[y][x]=min(dis.a[x][y],e[i].v);
}
st=lower_bound(p+1,p+n+1,st)-p,ed=lower_bound(p+1,p+n+1,ed)-p;
dis=quick_power(dis,k);
printf("%d\n",dis.a[st][ed]);
return 0;
}
浙公网安备 33010602011771号