「学习笔记」矩阵快速幂

一.快速幂#

P1226 【模板】快速幂||取余运算#

快速幂用来解决这样的式子 $a^b$，大家应该都会，我就放个代码。

Copy
ll ksm (ll a, ll b) {//a^b
	ll res = 1;
	while (b) {
		if (b & 1) {
			res = res * a;
		}
		b >>= 1;
		a = a * a;
	}
	return res;
}

Ybtoj【例题1】序列的第k个数#

这道题就可以直接根据前三项判断是等差数列还是等比数列。

等差数列 $a_n=a_1+(n-1)×q$，等比数列 $a_n=a_1×q_{n-1}$，$q$ 为公差。

---

二.矩阵#

概念1.单位矩阵#

主对角线上的元素都为 $1$，其余元素都为 $0$ 的 $n$ 阶矩阵称为 $n$ 阶单位矩阵。

单位矩阵的重要点在于：单位矩阵的任意次方都等于它本身。

如图，这就是一个三阶单位矩阵：

在矩阵快速幂转移时，可以将其设为初始矩阵。

---

概念2.矩阵加/减法#

两个矩阵相加/减，就将其对应位置上的数相加/减。

设 $A,B$ 为长宽相等的矩阵，矩阵 $C=A±B$，那么 $A_{i,j}±B_{i,j}=C_{i,j}$。

下图就是一个加法的例子。

---

概念3.矩阵乘法#

设矩阵 $C=A×B$，有几个要求。

$1.$ $A$ 的列数等于 $B$ 的行数。

$2.$ $A$ 是 $n×r$ 的矩阵， $B$ 是 $r×m$ 的矩阵，那么 $C$ 就是 $n×m$ 的矩阵。

$3.$ $C_{i,j}=\sum ^{r}_{k=1}A_{i,k}B_{k,j}$。

如图就是矩阵乘法的例子：

我们显然可以得到一个结论：

矩阵乘法满足结合律，不满足交换律。

这里给出 $O(n^3)$ 的矩阵乘法代码：

Copy
inline void Mul (int f[][], int x[][], int y[][]) {//计算矩阵x乘y=f。
	int t[N][N];//临时记录答案矩阵。
	for (int i = 0; i < n; i ++) {
		for (int j = 0; j < n; j ++) {
			t[i][j] = 0;//清空答案矩阵。
		}
	}
	for (int i = 0; i < n; i ++) {
		for (int j = 0; j < n; j ++) {
			for (int k = 0; k < n; k ++) {
				t[i][j] = t[i][j] + x[i][k] * y[k][j];//矩阵乘法法则。
			}
		}
	}
	for (int i = 0; i < n; i ++) {
		for (int j = 0; j < n; j ++) {
			f[i][j] = t[i][j];//复制答案。
		}
	}
}

---

概念4.方阵乘幂#

方阵是一个长和宽相等的矩阵。

将方阵 $A$ 自乘 $n$ 次，记作 $A^n$。

上文说到矩阵乘法满足结合律，那么我们就可以将方阵乘幂利用快速幂的思想优化。

---

三.矩阵快速幂的引入#

P3390 【模板】矩阵快速幂#

先解决这道模板题。

给定矩阵 $A$，求 $A^k$。

只需要将快速幂中的数字变为矩阵即可，乘法变为矩阵乘法。

给出一个模板，这是以后做题的基础。

Copy
inline void qpow (int T) {
    initans ();//这里是将ans初始化为单位矩阵。
    
    while (T) {
        if (T & 1) {
            Mul (ans, ans, a);
            //ans=ans*a。
        }
        
        Mul (a, a, a);
        //a=a*a。
        
        T >>= 1;
    }
}

观察一下，是不是与原来的快速幂模板一模一样？

---

Ybtoj【例题2】斐波那契数列#

P1962 斐波那契数列#

求斐波那契数列的第 $n$ 项 $\bmod10^9+7，n<2^{63}$。

$O(n)$ 递推显然不可做，我们就需要采用矩阵快速幂来进行优化。

我们的思路在于：

我们将最后求的答案变为一个矩阵 $B$，构建一个初始矩阵 $A$ 和一个转移矩阵 $base$，使 $A×base=B$。

那么，我们设答案矩阵 $fib(n)=\begin{bmatrix} f_{n} & f_{n-1} \end{bmatrix}$，我们希望通过 $fib(n-1)=\begin{bmatrix} f_{n-1} & f_{n-2} \end{bmatrix}$ 得到 $fib(n)$，也就是说让 $fib(n)=fib(n-1)×base$。

由 $f_n=f_{n-1}+f_{n-2}$ 得：

转移矩阵的第一列为 $\begin{bmatrix} 1 \\ 1 \end{bmatrix}$，因为 $f_{n-1}×1+f_{n-2}×1=f_n$。

那么可以推出转移矩阵的第二列为 $\begin{bmatrix} 1 \\ 0 \end{bmatrix}$，因为 $f_{n-1}×1+f_{n-2}×0=f_{n-1}$。

那么转移矩阵就是 $\begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}$。

因为 $f_1=f_2=1$，那么设初始矩阵 $fib(2)=\begin{bmatrix} 1 & 1 \end{bmatrix}$，乘一次 $base$ 都能计算一项，也就是说我们要从 $fib(2)$ 推到 $fib(n)$ 需要乘 $n-2$ 次 $base$。

也就说 $A×base^{n-2}=B$。

核心代码如下：

Copy
inline void qpow (int b) {
    while (b) {
        if (b & 1) {
            Mul (f, f, a);///计算f=f×a。
        }
        
        Mul (a, a, a);//a=a^2。
        
        b >>= 1;
    }
}

最后输出 $f_{0,0}$，答案矩阵的第一位就是 $f_n$。

这里的矩阵快速幂仅仅是把快速幂中的数字改为了矩阵而已。

---

P1939 【模板】矩阵加速（数列）#

设答案矩阵为 $F(n)=\begin{bmatrix} a_{n} & a_{n-1} & a_{n-2} \end{bmatrix}$，有一个矩阵 $F(n-1)=\begin{bmatrix} a_{n-1}&a_{n-1} & a_{n-3} \end{bmatrix}$，我们要构建矩阵 $base$ 使 $F(n-1)×base=F(n)$。

然后可以推出矩阵 $base=\begin{bmatrix} 1 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{bmatrix}$。

套模板即可。

---

四.矩阵快速幂的运用#

P1349 广义斐波那契数列#

易得转移矩阵 $base=\begin{bmatrix} p & 1 \\ q & 0 \end{bmatrix}$，初始矩阵 $A=\begin{bmatrix} a_{2} & a_{1} \end{bmatrix}$，套矩阵快速幂模板转移 $n-2$ 次即可。

P2044 [NOI2012] 随机数生成器#

题目给出了递推式，发现要加上 $c$，我们的初始矩阵就设为 $[x_0,c]$。

然后就根据递推式得出转移矩阵 $\begin{bmatrix} a & 0 \\ 1 & 1 \end{bmatrix}$，矩阵快速幂 $n$ 次转移即可。

但是我们发现数据非常大，那么就要使用龟速乘。

---

Ybtoj【例题3】行为方案#

P3758 [TJOI2017]可乐#

P5789 [TJOI2017]可乐（数据加强版）#

发现这个自爆操作不太好搞，可以创建一个虚拟点 $0$，自爆操作就相当于走到 $0$ 点。

可以设出一个状态 $f_{i,j,k}$ 为 $i$ 走到 $j$ 走 $k$ 步的方案数，那么就有一个暴力方程 $f_{i,j,k}=\sum$ $f_{i,p,k-1} × f_{p,j,1}$。

然后我们发现这个第三维可以压掉，那么 $f_{i,j}=\sum f_{i,p}×f_{p,j}$。

这样的复杂度是 $O(tn^2)$ 的，无法通过怎么办？

发现这个转移方程非常符合矩阵快速幂的形式，我们可以用矩阵快速幂优化到 $O(\log tn^2)$。

最后的答案就是 $\sum^{n}_{i=0} ans_{1,i}$。

代码如下：

Copy
#include <bits/stdc++.h>

using namespace std;

const int N = 105;
const int mod = 2017;

int n, m, t, base[N][N], ans[N][N];

inline void Mul (int f[N][N], int x[N][N], int y[N][N]) {
    int t[N][N];
    
    for (int i = 0; i <= n; i ++) {
        for (int j = 0; j <= n; j ++) {
            t[i][j] = 0;
        }
    } 
    
    for (int i = 0; i <= n; i ++) {
        for (int j = 0; j <= n; j ++) {
            for (int k = 0; k <= n; k ++) {
                t[i][j] = (t[i][j] + x[i][k] * y[k][j] % mod) % mod;
            }
        }
    }
    
    for (int i = 0; i <= n; i ++) {
        for (int j = 0; j <= n; j ++) {
            f[i][j] = t[i][j];
        }
    }
}

inline void qpow (int T) {
    while (T) {
        if (T & 1) {
            Mul (ans, ans, base);
        }
        
        Mul (base, base, base);
        
        T >>= 1;
    }
}

int main () {
    scanf ("%d%d", &n, &m);
    
    for (int i = 1; i <= m; i ++) {
        int u, v;
        
        scanf ("%d%d", &u, &v);
        
        base[u][v] = 1;
        base[v][u] = 1;//有连边的方案为1。
    }
    
    for (int i = 0; i <= n; i ++) {
        base[i][i] = 1;//自己到自己方案为1。
        base[i][0] = 1;//到爆炸点方案为1。
        ans[i][i] = 1;//初始矩阵设为单位矩阵。
    }
    
    scanf ("%d", &t);
    
    qpow (t);//因为初始的时刻为0，所以总共要转移t次。
    
    int sum = 0;
    
    for (int i = 0; i <= n; i ++) {
        sum = (sum + ans[1][i]) % mod;
    }
    
    printf ("%d", sum);
    
    return 0;
}

代码中初始矩阵和答案矩阵是同一个矩阵，因为我在快速幂时不断覆盖初始矩阵。

初始矩阵因为是没操作的情况，那么就设为单位矩阵。

---

Ybtoj【例题4】矩阵求和#

这题我就稍微讲一下思路。（因为我也不太会做）

我们构造矩阵 $B=\begin{bmatrix} A & I \\ 0 & I \end{bmatrix}$，其中 $A$ 为给出的初始矩阵， $I$ 为 $n×n$ 的单位矩阵， $0$ 为 $n×n$ 的零矩阵。

那么 $B^{k+1}=\begin{bmatrix} A^{k+1} & I+\Sigma _{i=1}^{k}A^{i} \\ 0 & 1 \end{bmatrix}$，就发现 $B^{k+1}$ 右上角的矩阵减去一个单位矩阵就是我们所要求的的答案。

---

Ybtoj【例题5】最短路径#

P2886 [USACO07NOV]Cow Relays G#

我们发现，有 $100$ 条边，但有 $1000$ 个点，那我们先离散化一下。

抛开边数限制，我们可以采用 $Floyd$ 算法。

图的邻接矩阵存储本质上就是一个矩阵，那我们可以想到用矩阵快速幂优化。

由于矩阵乘法满足结合律，我们直接转化，将乘法变为取 $\min$。

每次乘幂都是相当于走了一条边，那么乘 $n$ 次转移矩阵就相当于走了 $n$ 条边。

模板转化即可。

Copy
inline void Mul (int f[N][N], int x[N][N], int y[N][N]) {
    int t[N][N];
    for (int i = 0; i < tot; i ++) {
        for (int j = 0; j < tot; j ++) {
            t[i][j] = inf;//floyd算法初始为inf。
        }
    } 
    for (int i = 0; i < tot; i ++) {
        for (int j = 0; j < tot; j ++) {
            for (int k = 0; k < tot; k ++) {
                t[i][j] = min (t[i][j], x[i][k] + y[k][j]);//将乘法改为min。
            }
        }
    }
    for (int i = 0; i < tot; i ++) {
        for (int j = 0; j < tot; j ++) {
            f[i][j] = t[i][j];
        }
    }
}

完整代码。

多说一句，由于 $Ybtoj$ 上的数据有问题，所以有一个点需要数据特判才能 $AC$。

---

P2233 [HNOI2002]公交车路线#

这种题一看就是线性递推，那么就容易想到用矩阵快速幂进行优化。

阅读题面先写出几个信息：

$1.$ $n$ 为偶数时方案数为 $0$。

$2.$ $n<4$ 时方案数为 $0$。

那么我们先特判 $n$，然后把 $n$ 除以二，这样每一个数字都是相邻的，便于处理。

我们可以马上秒一个 $dfs$ 爆搜，得到前几位：
$f_1=0,f_2=2,f_3=8,f_4=28$。

猜测 $f_i$ 与 $f_{i-1}$ 和 $f_{i-2}$ 有关系，那我们列方程组，解得递推式如下：

$f_i=4×f_{i-1}-2×f_{i-2}$。

设答案矩阵为 $F(n)$，那么这样构建初始矩阵 $F(2)=[0, 2]$，转移矩阵 $base=\begin{bmatrix} 4 & 1 \\ -2 & 0 \end{bmatrix}$。

然后转移 $n-2$ 次就可以了。

答案要乘 $2$ ，因为我们只计算走向一边的方案。

又因为转移矩阵中有负数，所以我们要在答案输出时防止出现负数，像这样操作：
cout << (2 * ans[0][0] % mod + mod) % mod;

完整代码不放了就套模板即可。