矩阵快速幂（矩阵加速）

//南昌理工ACM集训队

本人小白如有不对欢迎指正ლ(╹◡╹ლ)

矩阵加速

快速幂（前置技能）

对于普通的求a的b次方，一般我们用心算pow(a , b)就可以立马得到想要的数据，可是当b的值为1e8、1e9的时候，普通的乘法就无能为力了，因为乘一次花费的时间太长了，大大增加了时间成本。所以我们如何解决这个问题呢？

我们知道，任何一个正整数都可以由2的整数次幂相加得到，换句话说就是每个数都可以转化为二进制（不说也知道吧），所以我们发现只要b的2进制的第i位为1，就乘上$a^{2^i}$ (这个值可以在循环中得到)

放个模板，如果不是很懂可以搜一搜别的博客，就不细讲了

#define ll long long
ll power(ll a, ll b)//a^b
{
    ll ans = 1;
    while (b)
    {
        if (b % 2 == 1) //当n为奇数时，乘以余下的一个a
            ans *= a;
        a *= a;
        b /= 2;
    }
    return ans;
}

矩阵快速幂

那么我们能否将矩阵代入快速幂的方法
有三个矩阵 a * b = c，
我们知道矩阵与矩阵相乘，其中c[i][j]为A的第i行与B的第j列对应乘积的和，即
$$cij=\sum _{k=1}^{n}aik\ast bkj$$

暴力循环小模板

juzhen mult(juzhen a, juzhen b) {
	juzhen c;
	for (int i = 0; i <= N; i++) {//N为矩阵的大小
		for (int j = 0; j <= N; j++) {
			for (int k = 0; k <= N; k++) {
				c.m[i][j] = (c.m[i][j] + a.m[i][k] * b.m[k][j]);
			}
		}
	}
	return c;
}

带进前面的快速幂的模板里，就是矩阵快速幂一开始的样子

juzhen result;
void juzhenPower(juzhen A, ll power) {
	while (power) {
		if (power & 1)result = mult(result, A);
		power /= 2;
		A = mult(A, A);
	}
}

当然，根据性质，1乘任何数都等于那个数，所以我们初始化的result矩阵乘任何矩阵都应该等于那个矩阵，所以根据矩阵的性质，我们可以用形如
$$\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)$$
的矩阵，就是左上到右下对角线都为1，其他数为0的矩阵。

那么矩阵快速幂该如何使用呢

可乐

洛谷传送门
acwing传送门（yxcyyds）

这题在洛谷也有数据减弱的版本可以用dp水过，当然我建议一步到位，直接矩阵加速ac两题，双倍快乐（

题意：有N个城市，给出M条城市之间的道路，有一个机器人有随机三种操作
① 停在原地
② 去下一个相邻的城市
③ 自爆
给出t时间，每秒机器人执行一种操作，问t秒内机器人有多少种操作。

对于 100% 的数据， n,m≤100 ， t≤10^9。

我们不难发现，n与m不是很大，但是t却大的离谱，建议机器人原地自爆

在大学的离散数学课中，我们学过用矩阵构造路径，去求有多少种路径的方法。在这题里，n与m最大只有100，也就是边长为100的矩阵，t为10 ^ 9,也就是将这个矩阵乘上10 ^ 9次方。利用矩阵快速幂的方法，我们也就可以求出我们需要的值。

那么我们按照构造路径的方法构造这个矩阵
设基础矩阵为base[ N ][ N ],base[ i ][ j ]表示从 i 到 j有多少种方法。当然最后求出的值就是base[ 1 ]上所有值的和

将样例

3 2
1 2
2 3
2

代入，可以构造出类似
$$\left(\begin{array}{cccc}1& 0& 0& 0\\ 1& 1& 1& 0\\ 1& 1& 1& 1\\ 1& 0& 1& 1\end{array}\right)$$
的矩阵，第0列的意思为从 i 到 0 ，也就是机器人的自爆，第 i 到 i 的意思为机器人站着不动，也就是花了一秒时间罚站。那么我们只要把这个矩阵乘 t 次方，就可以得到最后的记录方法数量的矩阵，将base[ 1 ]上所有的数相加，就能得到最终的答案

ac代码

#include<iostream>
#include<cstring>
using namespace std;
#define ll long long
const int mod = 2017;
ll N, M;
struct juzhen {
	ll m[101][101];
}x;
juzhen base, result;
juzhen az(juzhen a, juzhen b) {//暴力矩阵相乘
	juzhen c;
	memset(c.m, 0, sizeof(c.m));
	for (int i = 0; i <= N; i++) {
		for (int j = 0; j <= N; j++) {
			for (int h = 0; h <= N; h++) {
				c.m[i][j] = (c.m[i][j] + a.m[i][h] * b.m[h][j] % mod) % mod;
			}
		}
	}
	return c;
}
void juzhenPower(juzhen A, ll power) {//快速幂模板
	while (power) {
		if (power & 1)result = az(result, A);
		power /= 2;
		A = az(A, A);
	}
}
int main() {
	ll x, y, t, ans = 0;
	cin >> N >> M;
	for (int i = 0; i <= N; i++) {//将基础矩阵初始化
		base.m[i][i] = 1;
		base.m[i][0] = 1;
	}
	for (int i = 1; i <= N; i++) {//初始化result矩阵，也就是那个乘任意矩阵都等于那个矩阵的矩阵
		result.m[i][i] = 1;
	}
	for (int i = 0; i < M; i++) {//将连线代入到基础矩阵中
		cin >> x >> y;
		base.m[x][y] = 1;
		base.m[y][x] = 1;
	}
	cin >> t;
	juzhenPower(base, t);//将基础矩阵乘上t次
	for (int i = 0; i <= N; i++) {
		ans = (ans + result.m[1][i]) % mod;
	}
	cout << ans;//轻松ac
}

附个模板

矩阵相乘的部分有奆佬有时间复杂度为O( N^2.7)的方法，感兴趣的自己可以搜一搜

#define ll long long
juzhen cheng(juzhen a, juzhen b,ll mod) {//边长为N的两个矩阵相乘
	juzhen c;
	for (int i = 0; i <= N; i++) {
		for (int j = 0; j <= N; j++) {
			for (int h = 0; h <= N; h++) {
				c.m[i][j] = (c.m[i][j] + a.m[i][h] * b.m[h][j] % mod) % mod;
			}
		}
	}
	return c;
}
juzhen juzhenPower(juzhen A, ll n,ll mod) {//矩阵A的n次方，对mod取模，放心复制粘贴
	juzhen result;
	for (int i = 1; i <= N; i++) result.m[i][i] = 1;//初始化result
	while (n) {
		if (n & 1)result = cheng(result, A , mod);
		n /= 2;
		A = cheng(A, A, mod);
	}
	return result;
}

最后小结

当我们能将问题转化成求A矩阵的T次方时，我们就可以考虑能否用矩阵快速幂求出此题的答案，当然还要思考如何构造初始矩阵，当你能推出线性的状态转移方程时就可以根据状态转移方程构造初始矩阵，来解决类似的问题。