矩阵快速幂

8-10 upd：出大问题，矩阵顺序全写反了。晚点修。

9-27 upd：在修了。矩阵 LaTeX 好难写。

10-20 upd：修完了。

代码就先不管了。

---

前置知识：快速幂 & 矩阵乘法

快速幂

值得一提的是满足结合律的各种运算其实都可以用类似快速幂的这种每次对半折的方式优化。

矩阵乘法

分别给定 $n\times m$ 和 $m\times q$ 的矩阵 $A,B$，则两个矩阵相乘得到的 $n\times q$ 的矩阵 $C$ 满足：

$$C_{i,j}=\sum _{k=1}^m{A}_{i,k}\cdot B_{k,j}$$

举个例子，两个矩阵 $\left|\begin{array}{ccc}1& 2& 3\\ 4& 5& 6\end{array}\right|$，$\left|\begin{array}{cccc}7& 8& 9& 10\\ 11& 12& 13& 14\\ 15& 16& 17& 18\end{array}\right|$ 相乘得到的就是 $\left|\begin{array}{cccc}74& 80& 86& 92\\ 173& 188& 203& 218\end{array}\right|$，其中 $203=4\times 9+5\times 13+6\times 17$。

草率点说就是两个矩阵的行和列分别乘起来。用代码实现的话显然时间复杂度大概是立方级别的。

稍微再推几个例子就可以发现，矩阵乘法满足结合律。

但是要注意它并不满足交换律，所以在做有关题目的时候一定要注意计算顺序。

来自 8 月的 CT：你看这个垃圾博主就没注意顺序结果全文公式重构哈哈哈哈。

---

矩阵快速幂

实现

前面提到快速幂的思想还可用于其它满足结合律的运算，而矩阵乘法满足结合律。

所以当乘法中矩阵的行列数相等时，把这俩东西缝合一下就有了矩阵快速幂。非常简单，同样是对半折，多出折不掉的就给答案乘上，把普通快速幂里面的数字换成矩阵，把普通乘法换成矩阵乘法即可。设矩阵的边长为 $n$，幂的指数为 $k$，一次矩阵乘法 $O(n^3)$，所以一次矩阵快速幂就是 $O(n^3\log k)$ 的。

需要提到的是，普通快速幂里我们用于计算结果的变量最初会赋值为 $1$，而在矩阵快速幂里我们也显然需要一个这样的东西，不然如果拿个空的矩阵去乘显然啥都出不来。它就是单位矩阵。单位矩阵除了从左上角到右下角的对角线上全是 $1$ 之外其它位置都是 $0$，也就是 $\left|\begin{array}{cccc}1& 0& \cdots & 0\\ 0& 1& \cdots & 0\\ \cdots & \cdots & \cdots & \cdots \\ 0& 0& \cdots & 1\end{array}\right|$ 的形式，可以发现它乘任何矩阵都等于这个矩阵本身，类似普通乘法里的 $1$。

这里给出模板题的一种写法。

#include <bits/stdc++.h>
#define ll long long
#define pb push_back
#define pii pair<int,int>
using namespace std;
const ll N=107,MOD=1e9+7;
ll n,k;
struct mtr {
	ll m[N][N];
	mtr() { memset(m,0,sizeof(m)); } // 定义时清空数组
	void bui() { memset(m,0,sizeof(m)); for (int i=1;i<=n;i++) m[i][i]=1; } // 建出单位矩阵
} a,ans;
mtr operator * (mtr x,mtr y) { // 给矩乘重载个运算符
	mtr ret;
	for (ll i=1;i<=n;i++) for (ll j=1;j<=n;j++) for (ll l=1;l<=n;l++)
		ret.m[i][j]=(ret.m[i][j]+x.m[i][l]*y.m[l][j]%MOD+MOD)%MOD;
	return ret;
}
int main() {
	scanf("%lld%lld",&n,&k);
	for (ll i=1;i<=n;i++) for (ll j=1;j<=n;j++) scanf("%lld",&a.m[i][j]);
	ans.bui(); // 把这个存答案的矩阵开成单位矩阵
	while (k) {
		if (k&1) ans=ans*a;
		a=a*a,k>>=1;
	}
	for (ll i=1;i<=n;i++,printf("\n"))
		for (ll j=1;j<=n;printf("%lld ",ans.m[i][j]),j++);
	return 0;
}

应用

第一次学的人肯定会感觉矩阵幂这东西好像没啥用，但实际上用处可大了。

扔几个题。

P1939 矩阵加速

看这 $n$ 的范围，请出矩阵快速幂优化。

矩阵乘法的过程是行乘上列，那么我们如果往一个矩阵上摆一行变量，另一个矩阵上摆一列系数，那么一乘我们就可以得到一个新数。

于是根据这一点，$x\ge 4$ 时，根据 $a_x=a_{x-1}+a_{x-3}$ 这个递推式可以写出下面的矩阵乘法递推式：

$$\left|\begin{array}{ccc}{a}_x& a_{x-1}& a_{x-2}\end{array}\right|=\left|\begin{array}{ccc}{a}_{x-1}& a_{x-2}& a_{x-3}\end{array}\right|\times \left|\begin{array}{ccc}1& 1& 0\\ 0& 0& 1\\ 1& 0& 0\end{array}\right|$$

根据矩乘的过程理解一下。

构造这种矩阵的方法，草率点说就是把转移要用到的元素全部摆上去（比如这里要用到 $a_x,a_{x-1},a_{x-3}$，所以再补个 $a_{x-2}$ 全都放上去转移），$a$ 所在矩阵每一行对应转移矩阵的每一列，对应元素有多少就写多少（比如 $a_x=a_{x-1}+a_{x-3}=1\cdot a_{x-1}+0\cdot a_{x-2}+1\cdot a_{x-3}$，所以转移矩阵的第一列是 $101$，第二、三行同理。顺带一提到后面的题转移矩阵里放的就不一定全是 $0$ 和 $1$ 了）。

再把中间的 $\left|\begin{array}{ccc}{a}_{x-1}& a_{x-2}& a_{x-3}\end{array}\right|$ 层层代换，我们就可以得到：

$$\left|\begin{array}{ccc}{a}_x& a_{x-1}& a_{x-2}\end{array}\right|=\left|\begin{array}{ccc}{a}_3& a_2& a_1\end{array}\right|\times {\left|\begin{array}{ccc}1& 1& 0\\ 0& 0& 1\\ 1& 0& 0\end{array}\right|}^{x-3}=\left|\begin{array}{ccc}1& 1& 1\end{array}\right|\times {\left|\begin{array}{ccc}1& 1& 0\\ 0& 0& 1\\ 1& 0& 0\end{array}\right|}^{x-3}$$

所以直接用矩阵快速幂求出这个矩阵幂的值，再输出答案 $a_n$ 对应的结果即可。

代码的矩乘是反的，参考一下写法就差不多了。

#include <bits/stdc++.h>
#define ll long long
#define pb push_back
#define pii pair<int,int>
using namespace std;
const ll MOD=1e9+7;
ll n;
struct mtr {
	ll m[7][7];
	mtr() { memset(m,0,sizeof(m)); }
	void bui() { for (int i=1;i<=3;i++) m[i][i]=1; }
} a,ans;
mtr operator * (mtr r1,mtr r2) {
	mtr ret;
	for (ll i=1;i<=3;i++) for (ll j=1;j<=3;j++) for (ll l=1;l<=3;l++)
		ret.m[i][j]=(ret.m[i][j]+r1.m[i][l]*r2.m[l][j]%MOD)%MOD;
	return ret;
}
void solve() {
	scanf("%lld",&n);
	if (n<=3) printf("1\n");
	else {
		n-=3;
		memset(a.m,0,sizeof(a.m)),memset(ans.m,0,sizeof(ans.m)),ans.bui();
		a.m[1][1]=a.m[1][3]=a.m[2][1]=a.m[3][2]=1;
		while (n) {
			if (n&1) ans=ans*a;
			a=a*a,n>>=1;
		}
		printf("%lld\n",(ans.m[1][1]+ans.m[1][2]+ans.m[1][3])%MOD);
	}
}
int main() {
	int t; scanf("%d",&t); for (;t--;) solve();
	return 0;
}

P1349 广义斐波那契数列

同样构造矩阵，不过这次不只是 $0$ 和 $1$：

$$\left|\begin{array}{cc}{a}_n& a_{n-1}\end{array}\right|=\left|\begin{array}{cc}{a}_{n-1}& a_{n-2}\end{array}\right|\times \left|\begin{array}{cc}p& 1\\ q& 0\end{array}\right|$$

得到

$$\left|\begin{array}{cc}{a}_n& a_{n-1}\end{array}\right|=\left|\begin{array}{cc}{a}_2& a_1\end{array}\right|\times {\left|\begin{array}{cc}p& 1\\ q& 0\end{array}\right|}^{n-2}$$

同样输出 $a_n$ 对应的值即可。

#include <bits/stdc++.h>
#define ll long long
#define pb push_back
#define pii pair<int,int>
using namespace std;
ll p,q,a1,a2,n,MOD;
struct mtr {
	ll m[3][3];
	mtr() { memset(m,0,sizeof(m)); }
	void bui() { for (int i=1;i<=3;i++) m[i][i]=1; }
} a,ans;
mtr operator * (mtr r1,mtr r2) {
	mtr ret;
	for (ll i=1;i<=2;i++) for (ll j=1;j<=2;j++) for (ll l=1;l<=2;l++)
		ret.m[i][j]=(ret.m[i][j]+r1.m[i][l]*r2.m[l][j]%MOD)%MOD;
	return ret;
}
int main() {
	scanf("%lld%lld%lld%lld%lld%lld",&p,&q,&a1,&a2,&n,&MOD);
	if (n==1) printf("%lld",a1);
	else if (n==2) printf("%lld",a2);
	else {
		n-=2,ans.bui(),a.m[1][1]=p,a.m[2][1]=1,a.m[1][2]=q;
		while (n) {
			if (n&1) ans=ans*a;
			a=a*a,n>>=1;
		}
		printf("%lld",(ans.m[1][1]*a2%MOD+ans.m[1][2]*a1%MOD)%MOD);
	}
	return 0;
}

斐波那契数列前 $n$ 项和

题意：已知斐波那契数列 $f_1=f_2=1,f_3=2,f_4=3,\cdots ,f_n=f_{n-1}+f_{n-2}$。给定 $n,m$，求 $\sum _{i=1}^n{f}_i$ 对 $m$ 取模的值。$1\le n\le 2\times {10}^9$，$1\le m\le {10}^9+10$。

设 $S_i=\sum _{j=1}^i{f}_j$，则 $S_i=S_{i-1}+f_i=S_{i-1}+f_{i-1}+f_{i-2}=S_{i-1}+(S_{i-1}-S_{i-3})=2S_{i-1}-S_{i-3}$。

根据这个造出 $n\ge 4$ 时的矩阵递推式：

$$\left|\begin{array}{ccc}{S}_n& S_{n-1}& S_{n-2}\end{array}\right|=\left|\begin{array}{ccc}{S}_{n-1}& S_{n-2}& S_{n-3}\end{array}\right|\times \left|\begin{array}{ccc}2& 1& 0\\ 0& 0& 1\\ -1& 0& 0\end{array}\right|$$

得到

$$\left|\begin{array}{ccc}{S}_n& S_{n-1}& S_{n-2}\end{array}\right|=\left|\begin{array}{ccc}{S}_3& S_2& S_1\end{array}\right|\times {\left|\begin{array}{ccc}2& 1& 0\\ 0& 0& 1\\ -1& 0& 0\end{array}\right|}^{n-3}$$

当然还有一种做法是把 $f$ 和 $S$ 都扔到矩阵里转移：

$$\left|\begin{array}{ccc}{f}_n& f_{n-1}& S_n\end{array}\right|=\left|\begin{array}{ccc}{f}_{n-1}& f_{n-2}& S_{n-1}\end{array}\right|\times \left|\begin{array}{ccc}1& 1& 1\\ 1& 0& 1\\ 0& 0& 1\end{array}\right|$$

$$\left|\begin{array}{ccc}{f}_n& f_{n-1}& S_n\end{array}\right|=\left|\begin{array}{ccc}{f}_2& f_1& S_2\end{array}\right|\times {\left|\begin{array}{ccc}1& 1& 1\\ 1& 0& 1\\ 0& 0& 1\end{array}\right|}^{n-2}=\left|\begin{array}{ccc}1& 1& 2\end{array}\right|\times {\left|\begin{array}{ccc}1& 1& 1\\ 1& 0& 1\\ 0& 0& 1\end{array}\right|}^{n-2}$$

#include <bits/stdc++.h>
#define ll long long
#define pb push_back
#define pii pair<int,int>
#define pll pair<ll,ll>
#define vo void()
using namespace std;
ll n,m;
struct mtx {
	ll m[4][4];
	mtx() { memset(m,0,sizeof(m)); }
	void bui() { for (ll i=1;i<=3;i++) m[i][i]=1; }
} a,ans;
mtx operator * (mtx x,mtx y) {
	mtx ret;
	for (ll i=1;i<=3;i++) for (ll j=1;j<=3;j++) for (ll k=1;k<=3;k++)
		(ret.m[i][j]+=(x.m[i][k]*y.m[k][j]%m+m)+m)%=m;
	return ret;
}
int main() {
	scanf("%lld%lld",&n,&m);
	if (n==1) printf("1");
	else if (n==2) printf("2");
	else if (n==3) printf("4");
	else {
		ans.bui(),a.m[1][1]=2,a.m[2][1]=1,a.m[3][2]=1,a.m[1][3]=m-1,n-=3;
		while (n) {
			if (n&1) ans=ans*a;
			a=a*a,n>>=1;
		}
		printf("%lld",((ans.m[1][1]+m)*4+(ans.m[1][2]+m)*2+(ans.m[1][3]+m))%m);
	}
	return 0;
}

#include <bits/stdc++.h>
#define ll long long
#define pb push_back
#define pii pair<int,int>
#define pll pair<ll,ll>
#define vo void()
using namespace std;
ll n,m;
struct mtx {
	ll m[4][4];
	mtx() { memset(m,0,sizeof(m)); }
	void bui() { for (ll i=1;i<=3;i++) m[i][i]=1; }
} a,ans;
mtx operator * (mtx x,mtx y) {
	mtx ret;
	for (ll i=1;i<=3;i++) for (ll j=1;j<=3;j++) for (ll k=1;k<=3;k++)
		(ret.m[i][j]+=(x.m[i][k]*y.m[k][j]%m+m)+m)%=m;
	return ret;
}
int main() {
	scanf("%lld%lld",&n,&m);
	if (n==1) printf("1");
	else if (n==2) printf("2");
	else if (n==3) printf("4");
	else {
		ans.bui(),a.m[1][1]=a.m[1][2]=a.m[2][1]=a.m[3][1]=a.m[3][2]=a.m[3][3]=1,n-=2;
		while (n) {
			if (n&1) ans=ans*a;
			a=a*a,n>>=1;
		}
		printf("%lld",(ans.m[3][1]+ans.m[3][2]+ans.m[3][3]*2)%m);
	}
	return 0;
}

从这三道题我们不难看出，很多线性递推问题都是可以用矩阵快速幂优化到 $\log$ 级别的，优化方法也大同小异。

所以以后遇到线性递推问题我们就可以用矩阵快速幂薄纱它们了！

举个例子。

不知道什么东西

题意：$T$ 组数据。有 $N$ 块砖排成一排染色，每一块砖需要涂上红、蓝、绿、黄四种颜色中的其中一种，且当这 $N$ 块砖中红色和绿色的块数均为偶数时，染色效果最佳。问有多少种染色方案可以使效果达到最佳，答案对 $10007$ 取模。$1\le T\le 100$，$1\le N\le {10}^9$。

先不管数据范围，我们考虑用比较普通的方法怎么做。

设 $d{p}_{i,0}$ 为前 $i$ 块砖的染色方案中红绿砖块数量都为偶数的方案数，$d{p}_{i,1}$ 为红绿都为奇数的方案数，$d{p}_{i,2}$ 为红绿两种颜色有一种颜色的砖块数量是奇数的方案数。则不难得出 $d{p}_{i,0}=2\cdot d{p}_{i-1,0}+d{p}_{i-1,2}$，$d{p}_{i,1}=2\cdot d{p}_{i-1,1}+d{p}_{i-1,2}$，$d{p}_{i,2}=2\cdot d{p}_{i-1,0}+2\cdot d{p}_{i-1,1}+2\cdot d{p}_{i-1,2}$，其实就是把第 $i$ 块砖的所有颜色选择都砸进去。初始状态为 $d{p}_{i,0}=1$，$d{p}_{i,1}=d{p}_{i,2}=0$。

然后可以发现它是个线性递推，于是使用矩阵快速幂优化一下就可以过掉这题了：

$$\left|\begin{array}{ccc}d{p}_{n,0}& d{p}_{n,1}& d{p}_{n,2}\end{array}\right|=\left|\begin{array}{ccc}d{p}_{n-1,0}& d{p}_{n-1,1}& d{p}_{n-1,2}\end{array}\right|\times \left|\begin{array}{ccc}2& 0& 2\\ 0& 2& 2\\ 1& 1& 2\end{array}\right|$$

$$\left|\begin{array}{ccc}d{p}_{n,0}& d{p}_{n,1}& d{p}_{n,2}\end{array}\right|=\left|\begin{array}{ccc}d{p}_{0,0}& d{p}_{0,1}& d{p}_{0,2}\end{array}\right|\times {\left|\begin{array}{ccc}2& 0& 2\\ 0& 2& 2\\ 1& 1& 2\end{array}\right|}^n=\left|\begin{array}{ccc}1& 0& 0\end{array}\right|\times {\left|\begin{array}{ccc}2& 0& 2\\ 0& 2& 2\\ 1& 1& 2\end{array}\right|}^n$$

#include <bits/stdc++.h>
#define ll long long
#define pb push_back
#define pii pair<int,int>
#define pll pair<ll,ll>
#define vo void()
using namespace std;
const ll MOD=1e4+7;
struct mtx {
	ll m[4][4];
	mtx() { memset(m,0,sizeof(m)); }
	void bui() { memset(m,0,sizeof(m)); for (ll i=1;i<=3;i++) m[i][i]=1; }
} a,ans;
mtx operator * (mtx x,mtx y) {
	mtx ret;
	for (ll i=1;i<=3;i++) for (ll j=1;j<=3;j++) for (ll k=1;k<=3;k++)
		(ret.m[i][j]+=x.m[i][k]*y.m[k][j]%MOD)%=MOD;
	return ret;
}
void solve() {
	ll n;
	scanf("%lld",&n);
	ans.bui(),memset(a.m,0,sizeof(a.m));
	a.m[1][1]=2,a.m[1][3]=1,a.m[2][2]=2,a.m[2][3]=1,a.m[3][1]=a.m[3][2]=a.m[3][3]=2;
	while (n) {
		if (n&1) ans=ans*a;
		a=a*a,n>>=1;
	}
	printf("%lld\n",ans.m[1][1]);
}
int main() {
	int t; for (scanf("%d",&t);t--;solve());
	return 0;
}

接下来上点不怎么像线性递推的东西。

又是不知道什么东西

题意：多组数据。给定一张 $n$ 个点 $m$ 条边，点从 $0$ 到 $(n-1)$ 编号的有向图。有 $T$ 组询问，每次询问给定两个点 $s,t$，问存在多少条以 $s$ 为起点，$t$ 为终点恰好经过 $k$ 条边（可以重复）的不同路径。答案对 ${10}^3$ 取模。$1\le T\le 100$，$1\le n\le 20$，$0\le m\le 100$，$0\le k<20$。

答案可以很大，显然不能直接搜索。

非常神奇。

根据给定的边造出这张图的邻接矩阵，然后这个矩阵的 $k$ 次幂就是每对点之间恰好经过 $k$ 条边的路径数量，根据询问输出即可。

为什么？

回顾一下矩阵乘法的过程，可以发现这个乘法跟 Floyd 的松弛操作特别像。每乘一次，就相当于 Floyd 的把所有点对松弛一遍。只不过 Floyd 一般是求最短路，而这里是计数，所以是乘起来。

#include <bits/stdc++.h>
#define ll long long
#define pb push_back
#define pii pair<int,int>
using namespace std;
const ll N=27,MOD=1e3;
ll n,m,t;
struct mtr {
	ll m[N][N];
	mtr() { memset(m,0,sizeof(m)); }
	void bui() { memset(m,0,sizeof(m)); for (int i=0;i<n;i++) m[i][i]=1; }
} a,b,ans;
mtr operator * (mtr x,mtr y) {
	mtr ret;
	for (ll i=0;i<n;i++) for (ll j=0;j<n;j++) for (ll l=0;l<n;l++)
		ret.m[i][j]=(ret.m[i][j]+x.m[i][l]*y.m[l][j]%MOD)%MOD;
	return ret;
}
void solve() {
	memset(a.m,0,sizeof(a.m));
	for (ll i=1,x,y;i<=m;i++) scanf("%lld%lld",&x,&y),a.m[x][y]=1;
	scanf("%lld",&t);
	for (ll x,y,z;t--;) {
		scanf("%lld%lld%lld",&x,&y,&z);
		ans.bui(),b=a;
		while (z) {
			if (z&1) ans=ans*b;
			b=b*b,z>>=1;
		}
		printf("%lld\n",ans.m[x][y]);
	}
}
int main() {
	for (;scanf("%lld%lld",&n,&m),n+m;solve());
	return 0;
}

P2886 [USACO07NOV]Cow Relays

路径计数换成最短路照样可以用矩阵快速幂，只是这次的幂跟我们平常的认知不太一样……

之前写过，满足结合律的运算可以用快速幂优化。

实际上，不难发现，Floyd 求最短路的松弛过程 f[i][j]=min(f[i][j],f[i][k]+f[k][j]) 同样满足结合律！

所以我们把快速幂里的矩阵乘法改成 Floyd 的矩阵松弛就可以了，直接对乘法实现动刀。

不过这题比较恶心，点的编号稍微大了点（$1000$），但是边数只有 $100$，最终连得到边的点最多也就 $101$ 个（因为是连通图），所以只给连到了边的点加个离散化跑矩阵快速“幂”即可。

#include <bits/stdc++.h>
#define ll long long
#define pb push_back
#define pii pair<int,int>
#define pll pair<ll,ll>
#define vo void()
using namespace std;
const int N=107,M=1007;
int n,t,s,e,vis[M],cnt[M],cntn;
struct mtr { int m[N][N]; } ans,a;
mtr operator * (mtr x,mtr y) {
	mtr ret;
	for (int i=1;i<=cntn;i++) for (int j=1;j<=cntn;j++) ret.m[i][j]=1e9+1e7;
	for (int l=1;l<=cntn;l++)
		for (int i=1;i<=cntn;i++) for (int j=1;j<=cntn;j++)
			ret.m[i][j]=min(ret.m[i][j],x.m[i][l]+y.m[l][j]);
	return ret;
}
int main() {
	scanf("%d%d%d%d",&n,&t,&s,&e);
	for (int i=1;i<N;i++) for (int j=1;j<N;j++)
		ans.m[i][j]=a.m[i][j]=1e9+1e7;
	for (int i=1;i<N;i++) ans.m[i][i]=0;
	for (int i=1,x,y,z;i<=t;i++) {
		scanf("%d%d%d",&z,&x,&y);
		if (!vis[x]) cnt[x]=++cntn,vis[x]=1;
		if (!vis[y]) cnt[y]=++cntn,vis[y]=1;
		a.m[cnt[x]][cnt[y]]=a.m[cnt[y]][cnt[x]]=z;
	}
	while (n) { if (n&1) ans=ans*a; a=a*a,n>>=1; }
	printf("%d",ans.m[cnt[s]][cnt[e]]);
	return 0;
}

---

一点点拓展题

LOJ10222 佳佳的 Fibonacci：长得奇怪了点的递推优化。

AcWing225 矩阵幂的和：更奇怪的递推优化。

P4159 迷路：图上路径计数+奇妙建图小技巧。

P2151 HH去散步：图上路径计数+奇妙建图小技巧。

CF852B Neural Network Country：抽象的图上路径计数。

P3193 GT考试：利用 KMP 建图后路径计数。

题解啥的就先鸽了吧……