矩阵快速幂

矩阵快速幂是一种常用于DP的算法。通过矩阵乘法去快速转移状态求解。
但是由于矩阵运算的复杂度极高，因此一般而言矩阵快速幂优化DP的决策点不能太多。
矩阵快速幂的另一种应用是图论中的路径经过方案数的问题
具体来讲有一个结论：邻接矩阵的k次方表示某两个点之间路径距离为k的方案数，只不过需要注意在有向图中这种路径是有方向的。
当然，严格来讲这种矩阵快速幂仍然属于优化DP的范畴，因为其本质是对于floyd的扩展（当然扩展后也就和DP本身没什么关系了）
还有的题就是直接考察的矩阵本身，通过矩阵快速幂以及其它算法综合来解题。
接下来会具体分上面三种题型来总结矩阵快速幂相关的应用。

当然，在做矩阵快速幂优化前，先需要确保熟练掌握了矩阵快速幂本身P3390 【模板】矩阵快速幂

优化DP

P1962 斐波那契数列

斐波那契数列作为最经典的矩阵快速幂的题，一直作为其入门练习题的首选。
首先矩阵快速幂的两点性质一定要把握好：可DP，决策点固定且较少。而斐波那契数列就完美符合这两个要求。
然后就要考虑怎样去设计矩阵了。矩阵快速幂优化DP的主要思想是通过矩阵维护一些信息，使得可以通过这些信息可以推出下一个状态。
由于斐波那契数列是由上两个状态转移而来的，因此设矩阵

[\begin{matrix} f_{i} & f_{i - 1} \end{matrix}]

而下一个状态就可以通过这个矩阵转移。由于 $f_{i + 1} = f_{i} + f_{i - 1}$ ，因此只需要乘上一个简洁的矩阵就可以转移了。
这个过程稍显抽象，但看一下就懂了 ~~只可意会不可言传~~ 。对于 $f_{i + 1}$ ，有转移

[\begin{matrix} f_{i + 1} & f_{i} \end{matrix}] = [\begin{matrix} f_{i} & f_{i - 1} \end{matrix}] \times [\begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix}]

因此我们就构造出了转移矩阵以及具体的转移方法。直接快速幂就行了。
注意矩阵的码风很重要，一定要形成自己固定的写法。

一点行都没压

点击查看代码

#include<bits/stdc++.h>
using namespace std;
typedef unsigned long long ll;
#define int ll
const ll p=1e9+7;
ll n;
struct ju
{
	int a[3][3];
	ju(){memset(a,0,sizeof(a));}
}base,f;
ju operator * (const ju &x,const ju &y)
{
	ju tmp;
	for(int i=1;i<=2;i++)
		for(int j=1;j<=2;j++)
			for(int h=1;h<=2;h++)
				tmp.a[i][j]=(tmp.a[i][j]+x.a[i][h]*y.a[h][j]%p)%p;
	return tmp;
}
ju ksm(ju x,int k)
{
	ju res;
	for(int i=1;i<=2;i++) res.a[i][i]=1;
	while(k)
	{
		if(k&1) res=res*x;
		x=x*x;k>>=1;
	}
	return res;
}
signed main()
{
	cin>>n;
	if(n==1ll||n==2ll){cout<<'1';return 0;}
	base.a[1][1]=base.a[1][2]=base.a[2][1]=1,base.a[2][2]=0;
	ju ans=ksm(base,n-2);
	f.a[1][1]=f.a[1][2]=1;
	f=f*ans;
	cout<<f.a[1][1]<<'\n';
	return 0;
}

P1939 矩阵加速（数列）

此题本质上与上一道题是一样的，但是转移的位置变化了，递推式变为了 $f_{i} = f_{i - 1} + f_{i - 3}$
然后就来考虑一下转移的矩阵有什么变化。因为要从前一和三个位置转移，能不能直接维护这两个值呢？
事实上是不行的。自己手推一下就会发现，转移的时候，第一次转移确实能转，但之后矩阵里新的 $f_{i - 3}$ 就算不出来了。
所以考虑多维护一个东西。设答案矩阵为

[\begin{matrix} f_{i} & f_{i - 1} & f_{i - 2} \end{matrix}]

转移就有

[\begin{matrix} f_{i + 1} & f_{i} & f_{i - 1} \end{matrix}] = [\begin{matrix} f_{i} & f_{i - 1} & f_{i - 2} \end{matrix}] \times [\begin{matrix} 1 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{matrix}]

然后就与上面一模一样了。

点击查看代码

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=4;
const int p=1e9+7;
struct node
{
	int a[N][N];
	node(){memset(a,0,sizeof(a));}
}base,f;
node operator * (const node &x,const node &y)
{
	node res;
	for(int i=1;i<=3;i++)
	for(int j=1;j<=3;j++)
	for(int h=1;h<=3;h++) res.a[i][j]=(res.a[i][j]+x.a[i][h]*y.a[h][j])%p;
	return res;
}
node ksm(node x,int k)
{
	node res;for(int i=1;i<=3;i++) res.a[i][i]=1ll;
	while(k)
	{
		if(k&1) res=res*x;
		x=x*x,k>>=1;
	}
	return res;
}
signed main()
{
	ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
	int T;cin>>T;
	base.a[1][1]=base.a[1][2]=1;
	base.a[2][3]=1;
	base.a[3][1]=1;
	f.a[1][1]=f.a[1][2]=f.a[1][3]=1;
	while(T--)
	{
		int n;cin>>n;
		if(n<=3) {cout<<"1\n";continue;}
		node ans=ksm(base,n-3);
		ans=f*ans;
		cout<<ans.a[1][1]<<'\n';
	}
	return 0;
}

P3216 [HNOI2011] 数学作业

这道题就不是上面纯粹的递推了。

状态当然就是 $f_{i}$ 表示 $Concatenate (i) mod m$ 的值。现在来考虑一下转移方程。
如果学过对数，转移方程其实应该是比较好想的。有转移

f_{i} = (f_{i - 1} \times 10^{1 + \lg i} \mod m + i) \mod m

显然 $f_{i}$ 只与其前一个值有关，因此可以考虑矩阵快速幂去求解。答案与转移矩阵是显然的，有转移

[\begin{matrix} f_{i} & i + 1 & 1 \end{matrix}] = [\begin{matrix} f_{i - 1} & i & 1 \end{matrix}] \times [\begin{matrix} 10^{1 + \lg i} & 0 & 0 \\ 1 & 1 & 0 \\ 0 & 1 & 1 \end{matrix}]

由于 $\lg i$ 的总数很少，因此我们完全可以直接硬去枚举 $\lg i$ 的取值，复杂度仍然是正确的。最后将所有快速幂后的矩阵全部乘起来就可以了。
最后注意一下判断是否已经抵达了上界。

点击查看代码

#include<bits/stdc++.h>
using namespace std;
#define int unsigned long long
const int N=4;
int n,p;
struct node
{
	int a[N][N];
	node(){memset(a,0,sizeof(a));}
}base;
int lowbit(int x){return (-x)&x;}
node operator * (const node &x,const node &y)
{
	node res;
	for(int i=1;i<=3;i++)
		for(int j=1;j<=3;j++)
			for(int h=1;h<=3;h++)
			res.a[i][j]=(res.a[i][j]+x.a[i][h]*y.a[h][j]%p)%p;
	return res;
}
node ksm_ju(node x,int k)
{
	node ans=x;k--;
	while(k>0)
	{
		if(k&1) ans=ans*x;
		x=x*x;k>>=1;
	}
	return ans;
}
signed main()
{
	ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
	cin>>n>>p;
	base.a[2][1]=base.a[2][2]=base.a[3][1]=base.a[3][2]=base.a[3][3]=1;
	node ans;
	ans.a[1][3]=1;//ans.a[1][1]=1;
	for(int k=10;;k*=10)
	{
		base.a[1][1]=k%p;
		if(n<k) {ans=ans*ksm_ju(base,n-k/10+1);break;}
		ans=ans*ksm_ju(base,k-k/10);
	}
	cout<<ans.a[1][1];
	return 0;
}

图论中的应用

开头已经说过了，矩阵快速幂在图论中事实上与优化DP是类似的。这里先只做简要阐述。

P2886 [USACO07NOV] Cow Relays G

仍然是在开头，我们已经提到了一个结论：邻接矩阵的k次幂表示每一对点其经过k条边的路径种数。
我们考虑将这个结论类比一下用到这道题中来。
让我们先来思考一下上面的k次幂到底是怎样变成经过k条边的路径条数的？
发现比如对于1到2，如果k等于2，矩阵乘法时就会使比如1到3然后3到2的路径统计出来一一对应。
现在再来思考对于最短路问题，我们是否可以类似地去思考。事实上每做一次floyd，都是将两条最短路连起来，与上面类似。
同时，最短路本身也有类似于结合律的性质。分开求最短路以及拆分开求最短路是等价的。
所以我们改变一下矩阵乘法，将其变成求最短路，其代码与floyd很类似。当然，这里的求最短路只是简单地求min，与DP本身有些本质区别。

比较详细的就看代码吧

点击查看代码

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=2005;
int n,t,s,e,id[N],idcnt=0;
struct node
{
	int a[105][105];
	node(){memset(a,0,sizeof(a));}
}tu;
void putju(node x)
{
	for(int i=1;i<=idcnt;i++)
	{
		for(int j=1;j<=idcnt;j++)cout<<x.a[i][j]<<' ';
		cout<<'\n'; 
	}
}
node operator * (const node &x,const node &y)
{
	node res;memset(res.a,0x3f3f3f,sizeof(res.a));
	for(int i=1;i<=idcnt;i++)
		for(int j=1;j<=idcnt;j++)
			for(int h=1;h<=idcnt;h++) 
			res.a[i][j]=min(res.a[i][j],x.a[i][h]+y.a[h][j]);
	return res;
}
node ksm(node x,int k)
{
	node res;
	memset(res.a,0x3f3f3f,sizeof(res.a));
	for(int i=1;i<=idcnt;i++) res.a[i][i]=0;
	while(k)
	{
		if(k&1) res=res*x;
		x=x*x;k>>=1;
	}
	return res;
}
signed main()
{
	ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
	cin>>n>>t>>s>>e;
	memset(id,0,sizeof(id));memset(tu.a,0x3f3f3f,sizeof(tu.a));
	for(int i=1,u,v,w;i<=t;i++)
	{
		cin>>w>>u>>v;
		if(id[u]) u=id[u];else id[u]=++idcnt,u=id[u];
		if(id[v]) v=id[v];else id[v]=++idcnt,v=id[v];
		tu.a[u][v]=tu.a[v][u]=min(tu.a[u][v],w);
	}
	node ans=ksm(tu,n);
	cout<<ans.a[id[s]][id[e]]<<'\n';
	return 0;
}

其它类型

Power of Matrix

题意：多测，最多20组询问。给定矩阵 $A$ 求 $A + A^{2} + A^{3} \dots A^{k}, k \leq 1 e 6, n \leq 40$ ，其中n代表矩阵的边长。

显然不能直接去硬快速幂算。
发现一个比较显然的性质：如果我们已经求出来了 $A + A^{2} + \dots A^{\frac{k}{2}}$ ，就可以通过全部乘 $A^{\frac{k}{2}}$ 算出整体（当然k为奇数需要特判一下）
然后就发现这样可以一直分下去，我们实际上需要计算的只有 $l o g k$ 次矩阵乘法以及矩阵快速幂，因此就直接这样做就做完了，注意判边界。

点击查看代码

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=50;
const int p=10;
int n,m;
struct node
{
	int a[N][N];
	node() {memset(a,0,sizeof(a));}
}a;
node operator * (const node &x,const node &y)
{
	node res;
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
			for(int h=1;h<=n;h++)
			res.a[i][j]=(res.a[i][j]+x.a[i][h]*y.a[h][j]%p)%p;
	return res;
}
node operator + (const node &x,const node &y)
{
	node res;
	for(int i=1;i<=n;i++)
	for(int j=1;j<=n;j++) 
	res.a[i][j]=(x.a[i][j]+y.a[i][j])%p;
	return res;
}
node ksm(node x,int k)
{
	node res;for(int i=1;i<=n;i++) res.a[i][i]=1;
	while(k)
	{
		if(k&1) res=res*x;
		x=x*x,k>>=1;
	}
	return res;
}
node solve(node x,int k)
{
	if(k==1) return x;
	node tmp1=solve(x,k/2);
	node tmp2=tmp1*ksm(x,k/2);
	if(k%2==0) return tmp1+tmp2;
	else return tmp1+tmp2+ksm(x,k);
}
signed main()
{
	ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
	node ans;
	while(cin>>n>>m&&n)
	{
		memset(ans.a,0,sizeof(ans.a));
		memset(a.a,0,sizeof(a.a));
		for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++) cin>>a.a[i][j],a.a[i][j]%=p;
		ans=solve(a,m);
		for(int i=1;i<=n;i++)
		{
			for(int j=1;j<=n;j++) 
			if(j!=n)cout<<ans.a[i][j]<<' ';
			else cout<<ans.a[i][j];
			cout<<'\n';
		}
		cout<<'\n';
	}
	return 0;
}

P6569 [NOI Online #3 提高组] 魔法值

乍一看很不可做，但看一眼数据范围就知道是矩阵了。一共100个点，考虑如何构造矩阵去快速求解。
由于点数很少我们当然用邻接表去存，观察一下邻接表以及题目给出的这个式子，发现：对于每一个点，其需要按题目要求异或的所有值都是其邻接表里存的点。
同时手动模拟一下，就可以发现仍然对于邻接表的次幂，对于新算出来的这个矩阵的每一个位置就是相应需要异或的次数。
但乘法以及异或并没有结合律，因此我们不能直接去快速幂。
由于异或两次就抵消了，因此我们完全没有必要去做矩阵乘法。考虑更改一下矩阵乘法，将其“乘”变为

node res;res.h=x.h,res.l=y.l;
for(int i=1;i<=res.h;i++)
for(int j=1;j<=res.l;j++)
for(int k=1;k<=x.l;k++) res.mp[i][j]^=x.mp[i][k]*y.mp[k][j];

这样我们“乘”出来的矩阵每一位上就只有0或者1（因为本来就是邻接矩阵去“乘”，其初始时就只有01，那经过各种异或后仍然只有01）

显然对于一个01矩阵，某两个矩阵先乘起来后异或上这个01矩阵以及两个分别异或上这个01矩阵再乘起来是等价的。
因此现前给定的序列是 $f$ ，邻接矩阵是 $m p$ ，最终答案就是 $a n s_{a_{i}} = f x o r m p^{a_{i}}$ ，这里的 $x o r$ 就是指的上面的“乘”运算
但是这样复杂度是 $O (q n^{3} \sum_{i = 1}^{q} l o g a_{i})$ 的。发现主要影响复杂度的就是较多的矩阵乘法。 $O (n^{3})$ 实在是难以接受。
那如何去减少矩阵乘法的复杂度呢？发现每次矩阵乘法都会跑满 $O (n^{3})$ 的原因是我们是邻接表与邻接表相乘。
而我们发现答案矩阵行数只有1，那如果答案矩阵与邻接矩阵相乘那复杂度只有 $O (n^{2})$
因此考虑预处理出所有 $m p^{2^{k}}$ ，然后在回答询问时二进制拆位，这样就可以做到回答询问时的每次矩阵乘法只有 $O (n^{2})$ ，同时仍然只带一个 $l o g a_{i}$
这时总复杂度便达到了 $O (n^{3} l o g a_{m a x} + q n^{2} \sum_{i = 1}^{q} l o g a_{i})$ ，q、n同阶，就是 $O (n^{3} l o g a)$ 。（感觉比上一道难多了，最后的优化确实没想到）

事实上最后已经与快速幂没什么关系了，但是仍然是矩阵快速幂的思想。
同时其没有放在图论或者DP是因为其本身并没有明显的转移方程或者对图论性质的应用，而只是对矩阵本身优良性质、转化以及优化的考察。

点击查看代码

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=105;
const int M=(1ll<<33);
int n,m,q;
struct node
{
	int h,l,mp[N][N];
	node() {h=0,l=0,memset(mp,0,sizeof(mp));}
}f,yu[N];
node operator ^ (const node &x,const node &y)
{
	node res;res.h=x.h,res.l=y.l;
	for(int i=1;i<=res.h;i++)
	for(int j=1;j<=res.l;j++)
	for(int k=1;k<=x.l;k++) res.mp[i][j]^=x.mp[i][k]*y.mp[k][j];
	return res;
}
signed main()
{
	ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
	cin>>n>>m>>q;
	for(int i=1;i<=n;i++) cin>>f.mp[1][i];f.h=1,f.l=n;
	yu[0].h=yu[0].l=n;for(int i=1,u,v;i<=m;i++)cin>>u>>v,yu[0].mp[u][v]=yu[0].mp[v][u]=1;
	for(int i=1;i<=32;i++) yu[i].h=yu[i].l=n,yu[i]=yu[i-1]^yu[i-1];
	for(int i=1,x;i<=q;i++)
	{
		cin>>x;node res=f;res.h=1,res.l=n;
		for(int i=0;i<32;i++) if(x&(1<<i)) res=res^yu[i];
		cout<<res.mp[1][1]<<'\n';
	}
	return 0;
}

posted @ 2024-12-06 17:04 all_for_god 阅读(16) 评论(0) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

相关博文：

· 决策单调性DP

· 根号算法（分块、莫队、根号分治）

· 「学习笔记」矩阵快速幂

· 矩阵快速幂

· 数学/数论专题-学习笔记：矩阵小记#2（矩阵快速幂）

阅读排行：
· 分享一个免费、快速、无限量使用的满血 DeepSeek R1 模型，支持深度思考和联网搜索！
· 基于 Docker 搭建 FRP 内网穿透开源项目（很简单哒）
· ollama系列01：轻松3步本地部署deepseek，普通电脑可用
· 25岁的心里话
· 按钮权限的设计及实现

公告

昵称： all_for_god
园龄： 7个月
粉丝： 6
关注： 10

+加关注

2025年3月

日

一

二

三

四

五

六

allforgod

矩阵快速幂

优化DP

P1962 斐波那契数列

P1939 矩阵加速（数列）

P3216 [HNOI2011] 数学作业

图论中的应用

P2886 [USACO07NOV] Cow Relays G

其它类型

Power of Matrix

P6569 [NOI Online #3 提高组] 魔法值

公告

搜索

随笔分类

随笔档案

阅读排行榜

推荐排行榜

最新评论