文章目录
- 概述
- 快速幂
- 解析
- 代码
- 矩阵运算定义
- 加法
- 乘法
- 单位矩阵
- 一、斐波拉契（基础模板）
- 题目描述
- 解析
- 代码
- 二、行为方案（实际应用）
- 题目描述
- 解析
- 代码
- 三、矩阵求和（子矩阵作为矩阵元素）
- 题目描述
- 解析
- 代码
- 四、最短路径（对矩阵乘法灵活的定义方式）
- 题目描述
- 解析
- 代码：
- thanks for reading！
概述
矩阵快速幂，顾名思义，就是矩阵乘法与快速幂的结合,可以将时间复杂度优化到log级别。
当数据范围中图的尺寸较小，而时间、变换次数、步数等的范围较大（通常在10^9左右）时，可以考虑矩阵快速幂问题
矩阵快速幂的关键在于构造转移矩阵
快速幂
解析
由于
x^a * x^b = x^(a+b)
逆向来看,对于一个较大的指数k，我们可以进行拆分：
x^k = x^a1 * x^a2 *…(a1+a2+…=k)
由于二进制拆分的唯一存在性，我们就可以把k拆成若干个2的幂数作为a序列，将其对应的乘幂求出，并乘在一起就能得到x^k了
时间复杂度log k
代码
ll ksm(ll x, int w) {//求x的w次幂ll ans = 1;while (w) {if (w & 1)ans = (ll)(ans * x) % mod;//mod是取模用的x = (x * x) % mod;w >>= 1;}return ans;
}
矩阵运算定义
加法
对于两个大小完全相同的矩阵，定义二者相加为每个对应元素相加（比较简单，不细说了）
乘法
对于两个矩阵AB，若A的行数等于B的列数，则二者可以相乘
换言之，ab的矩阵A可以和bc的矩阵B相乘，结果是一个a*c的矩阵C，单次相乘时间复杂度为abc 的乘积
对于每一个c中的元素：
c[i][j]=∑a[i][p]+b[p][j](1<=p<=b)
简单说就是顶针式的乘积累加
容易证明，矩阵乘法满足结合律，但是不满足交换律
对于一个行列相等的方阵，可以不断累乘自身
单位矩阵
如果一个矩阵的结构为：
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1
则称它为一个单位矩阵
对于一个a*b的矩阵A，它与一个边长为b的矩阵相乘，结果还是A本身
单位矩阵在矩阵运算中起着单位1的作用
矩阵简单的运算就是这些，接下来我们通过一些例题由浅入深，看看矩阵快速幂如何应用吧
一、斐波拉契（基础模板）
万物之始
题目描述
斐波拉契的定义为：
f[1]=f[2]=1;
f[n]=f[n-1]+f[n-2](n>=3)
（不会真的有人不知道吧）
现在要你求出第k项取模后的结果
(k<=2^63)
解析
我们发现：
对于一个行矩阵：
f[n] f[n-1]
尝试构造出一个转移矩阵M：
（转移矩阵是矩阵题目的灵魂，多为方阵）
1 1
1 0
那么两者相乘后就会得到：
f[n]+f[n-1] f[n]
根据斐波拉契的定义，上面的结果就是：
f[n+1] f[n]
所以，要求第k项的值，就只需要把矩阵f[2] f[1]乘上k-2次转移矩阵即可
这样就可以用快速幂提速了
（由于本题原矩阵过于简单，其实代码实现时我们都不需要让那个行矩阵出现）
代码
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int mod=1e9+7;
ll a,b,c,k;
ll ans[4][4],res[4][4],trans[4][4];
void mul1(){//res*resmemset(trans,0,sizeof(trans));for(int i=1;i<=2;i++){for(int j=1;j<=2;j++){for(int p=1;p<=2;p++){trans[i][j]+=res[i][p]*res[p][j];trans[i][j]%=mod;}}}for(int i=1;i<=2;i++){for(int j=1;j<=2;j++){res[i][j]=trans[i][j];}}
}
void mul2(){//res*ansmemset(trans,0,sizeof(trans));for(int i=1;i<=2;i++){for(int j=1;j<=2;j++){for(int p=1;p<=2;p++){trans[i][j]+=ans[i][p]*res[p][j];trans[i][j]%=mod;}}}for(int i=1;i<=2;i++){for(int j=1;j<=2;j++){ans[i][j]=trans[i][j];}}
}
void ksm(ll w){while(w){if(w&1) mul2();mul1();w>>=1;}return;
}
void print(){for(int i=1;i<=2;i++){for(int j=1;j<=2;j++){printf("%d ",ans[i][j]);}printf("\n");}
}
int main()
{ans[1][1]=ans[2][2]=1;res[1][1]=res[1][2]=res[2][1]=1;scanf("%lld",&a);if(a<=2) {printf("1");return 0;}ksm(a-2);
// print();printf("%lld",(ans[1][1]+ans[2][1])%mod);return 0;
}二、行为方案（实际应用）
题目描述
解析
定义dp[i][j][t]为经过t时间从i走到j的方案数
不难想到转移
dp[i][j][t]=∑dp[i][p][t-1]*dp[p][j][1]
我们发现出现了矩阵乘法标志性的顶针结构
也就是说，想要从t-1转移到t，只需要把t-1状态的矩阵与状态1的矩阵相乘就行了
换言之，询问t状态的矩阵，就是求状态1矩阵的t次幂
这样在提速的同时，也就把第三维空间优化掉了
这样就转化为快速幂问题了
接下来就是原始矩阵的构造问题
首先肯定要把存在边的点连上：
dp[from][to]=1
还可以停留在原点，也就相当于一个自环：
dp[i][i]=1;(1<=i<=n)
对于自爆，我们可以把0结点当成自爆的状态，在任何一个点均可自爆，所以：
dp[i][0]=1(1<=i<=n)
由于自爆后不能向任何状态转移了，所以0结点没有任何出边
最后的答案就是∑dp[1][i] (0<=i<=n）
这样本题就结束了
代码
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int mod=2017;
const int N=150;
ll a,b,c,k;
int n,m;
ll ans[N][N],res[N][N],trans[N][N];
void mul1(){//res*resmemset(trans,0,sizeof(trans));for(int i=0;i<=n;i++){for(int j=0;j<=n;j++){for(int p=0;p<=n;p++){trans[i][j]+=res[i][p]*res[p][j];trans[i][j]%=mod;}}}for(int i=0;i<=n;i++){for(int j=0;j<=n;j++){res[i][j]=trans[i][j];}}
}
void mul2(){//res*ansmemset(trans,0,sizeof(trans));for(int i=0;i<=n;i++){for(int j=0;j<=n;j++){for(int p=0;p<=n;p++){trans[i][j]+=ans[i][p]*res[p][j];trans[i][j]%=mod;}}}for(int i=0;i<=n;i++){for(int j=0;j<=n;j++){ans[i][j]=trans[i][j];}}
}
void ksm(ll w){while(w){if(w&1) mul2();mul1();w>>=1;}return;
}
void print(){for(int i=0;i<=n;i++){for(int j=0;j<=n;j++){printf("%d ",ans[i][j]);}printf("\n");}
}
int main()
{scanf("%d%d",&n,&m);for(int i=0;i<=n;i++) ans[i][i]=1;for(int i=1;i<=m;i++){int a,b;scanf("%d%d",&a,&b);res[a][b]=res[b][a]=1;}for(int i=0;i<=n;i++) res[i][i]=1;for(int i=1;i<=n;i++) res[i][0]=1;int t;scanf("%d",&t);ksm(t);int tot=0;for(int i=0;i<=n;i++){tot+=ans[1][i];}
// print();printf("%d",tot%mod);return 0;
}三、矩阵求和（子矩阵作为矩阵元素）
题目描述
解析
对于矩阵A，我们构造转移矩阵：
| A | I |
|---|---|
| 0 | I |
（I表示与A边长相等的单位矩阵，矩阵套矩阵表示把小矩阵值套进大矩阵里）
比如，当A为：
2 3
1 5
时，转移矩阵为：
2 3 1 0
1 5 0 1
0 0 1 0
0 0 0 1
使它与自身相乘,得到：
| A^2 | I+A |
|---|---|
| 0 | I |
再乘一次：
| A^3 | I+A+A^2 |
|---|---|
| 0 | I |
这样我们就可以得到规律了
只需要把这个转移矩阵构造出来后乘上k次幂在把左上角减去单位矩阵输出即可
注意：取模减法需要判断一下会不会减成负的！
（连ybt评测的答案似乎都没有注意。。。）
代码
#include<bits/stdc++.h>
using namespace std;
#define ll long long
int mod=2017;
const int N=150;
ll a,b,c,k;
int n,m;
ll ans[N][N],res[N][N],trans[N][N];
void mul1(){//res*resmemset(trans,0,sizeof(trans));for(int i=1;i<=n;i++){for(int j=1;j<=n;j++){for(int p=1;p<=n;p++){trans[i][j]+=res[i][p]*res[p][j];trans[i][j]%=mod;}}}for(int i=1;i<=n;i++){for(int j=1;j<=n;j++){res[i][j]=trans[i][j];}}
}
void mul2(){//res*ansmemset(trans,0,sizeof(trans));for(int i=1;i<=n;i++){for(int j=1;j<=n;j++){for(int p=1;p<=n;p++){trans[i][j]+=ans[i][p]*res[p][j];trans[i][j]%=mod;}}}for(int i=1;i<=n;i++){for(int j=1;j<=n;j++){ans[i][j]=trans[i][j];}}
}
void ksm(ll w){while(w){if(w&1) mul2();mul1();w>>=1;}return;
}
void print(){for(int i=1;i<=n;i++){for(int j=1+n;j<=2*n;j++){if(i+n==j) ans[i][j]--;
// if(ans[i][j]<0) ans[i][j]+=mod;printf("%lld ",ans[i][j]);}printf("\n");}
}
int main()
{scanf("%d%d%d",&n,&m,&mod);for(int i=1;i<=n;i++){for(int j=1;j<=n;j++) scanf("%lld",&res[i][j]);}for(int j=n+1;j<=2*n;j++){res[j][j]=res[j-n][j]=1;}n<<=1;for(int i=1;i<=n;i++) ans[i][i]=1;ksm(m+1);n>>=1;print();return 0;
}四、最短路径（对矩阵乘法灵活的定义方式）
题目描述
解析
定义dp[i][j][k]：从i到j，经过k条边的最短路径
易得转移方程：
dp[i][j][k]=min(dp[i][j][k],dp[i][p][k-1]+dp[p][j][1];
我们发现出现了顶针结构，但似乎运算从乘积累加变成了求min
那么把运算的定义改一下不就好了！
（这样修改之后不再有单位矩阵的概念，好在本题也不需要）
还有一些细节问题：
1.在原始矩阵中dp[i][i]应该为正无穷而不是0，换言之，经过一条边走到原处的方案在没有自环的情况下应该是不可能的
2.由于边数<=100，所以点不会超过200个，但编号会达到1000。如果直接用原编号，会超时。所以本题需要对点进行离散化处理
代码：
#include<bits/stdc++.h>
using namespace std;
#define ll long long
int mod=2017;
const int N=1500;
int s,e,n,m,k;
int mx;
int a,b,c,d;
ll ans[N][N],res[N][N],trans[N][N];
int id[N*100],tot=0;
void mul1(){//res*resmemset(trans,0x3f,sizeof(trans));for(int i=1;i<=n;i++){for(int j=1;j<=n;j++){for(int p=1;p<=n;p++){trans[i][j]=min(trans[i][j],res[i][p]+res[p][j]);
// trans[i][j]%=mod;}}}for(int i=1;i<=n;i++){for(int j=1;j<=n;j++){res[i][j]=trans[i][j];}}
}
void mul2(){//res*ansmemset(trans,0x3f,sizeof(trans));for(int i=1;i<=n;i++){for(int j=1;j<=n;j++){for(int p=1;p<=n;p++){trans[i][j]=min(trans[i][j],ans[i][p]+res[p][j]);
// trans[i][j]%=mod;}}}for(int i=1;i<=n;i++){for(int j=1;j<=n;j++){ans[i][j]=trans[i][j];}}
}
void ksm(ll w){while(w){if(w&1) mul2();mul1();w>>=1;}return;
}
void print(){for(int i=1;i<=n;i++){for(int j=1;j<=n;j++){
// if(ans[i][j]<0) ans[i][j]+=mod;if(ans[i][j]>2e9) printf("-1 ");else printf("%lld ",ans[i][j]);}printf("\n");}printf("\n");
}
void print_res(){for(int i=1;i<=n;i++){for(int j=1;j<=n;j++){
// if(ans[i][j]<0) ans[i][j]+=mod;if(res[i][j]>2e9) printf("-1 ");else printf("%lld ",res[i][j]);}printf("\n");}
}
int main()
{scanf("%d%d%d%d",&k,&n,&s,&e);memset(res,0x3f,sizeof(res));memset(ans,0x3f,sizeof(ans));for(int i=1;i<=n;i++){scanf("%d%d%d",&c,&a,&b);if(!id[a]) id[a]=++tot;if(!id[b]) id[b]=++tot;a=id[a],b=id[b];res[a][b]=res[b][a]=c;
// mx=max(mx,max(a,b));
// printf("a=%d b=%d v=%d\n",a,b,c);}n=tot;
// print_res();for(int i=1;i<=n;i++) ans[i][i]=0;ksm(k);
// mul2();
// print();
// mul2();
// print();printf("%lld",ans[id[s]][id[e]]);return 0;
}
/*
1 6 6 4
11 4 6
4 4 8
8 4 9
6 6 8
2 6 9
3 8 9
*/
【插头dp】)
![牛客题霸 [将字符串转化为整数] C++题解/答案](http://pic.xiahunao.cn/牛客题霸 [将字符串转化为整数] C++题解/答案)
![P6800-[模板]Chirp Z-Transform【NTT】](http://pic.xiahunao.cn/P6800-[模板]Chirp Z-Transform【NTT】)
![牛客题霸 [滑动窗口的最大值] C++题解/答案](http://pic.xiahunao.cn/牛客题霸 [滑动窗口的最大值] C++题解/答案)
![牛客题霸 [在转动过的有序数组中寻找目标值] C++题解/答案](http://pic.xiahunao.cn/牛客题霸 [在转动过的有序数组中寻找目标值] C++题解/答案)
![牛客题霸 [栈和排序] C++题解/答案](http://pic.xiahunao.cn/牛客题霸 [栈和排序] C++题解/答案)
![邮局-[IOI2000]（四边形不等式）](http://pic.xiahunao.cn/邮局-[IOI2000]（四边形不等式）)
![P6847-[CEOI2019]Magic Tree【dp,线段树合并】](http://pic.xiahunao.cn/P6847-[CEOI2019]Magic Tree【dp,线段树合并】)