题解：P3390 【模板】矩阵快速幂 & 矩阵快速幂加速递推的应用思路

前置知识

今天把矩阵快速幂又复习了一遍，来博客上水一篇题解吧。
总而言之，矩阵快速幂 = 矩阵乘法 + 快速幂。

那矩阵乘法是什么呢？

设有两个矩阵 $ A $ 和 $ B $，其中 $ A $ 是一个 $ m \times n $ 矩阵，$ B $ 是一个 $ n \times p $ 矩阵。它们的乘积 $ C = A \times B $ 将是一个 $ m \times p $ 矩阵 $ C $，其中每个元素 $c_{ij} $ 的计算公式为：

\(c_{ij} = \sum_{k=1}^{n} a_{ik} b_{kj}\)

其中：

• $ c_{ij} $ 是矩阵 $ C $ 中第 $ i $ 行第 $ j $ 列的元素。
• $ a_{ik} $ 是矩阵 $ A $ 中第 $ i $ 行第 $ k $ 列的元素。
• $ b_{kj} $ 是矩阵 $ B $ 中第 $ k $ 行第 $ j $ 列的元素。
• $ n $ 是矩阵 $ A $ 的列数（也是矩阵 $ B $ 的行数）。

注意，我们有一点非常重要，我们能否把\(AB\)相乘取决于一个关键的条件:
\(A\)矩阵的列数必须等于\(B\)矩阵的行数。
因此：
\(AB \neq BA\)。

代码模板:

void mx(int a[N][N], int b[N][N], int n) {
    memset(tmp, 0, sizeof(tmp));
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            for (int k = 1; k <= n; k++) {
                tmp[i][j] += a[i][k] * b[k][j];
                tmp[i][j] %= m;
            }
        }
    }
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            a[i][j] = tmp[i][j];
        }
    }
}

快速幂：
用于快速计算 \(a^{b} mod m\)。（在 \(O(logn)\) 的复杂度下)。
oiwiki直链：https://oiwiki.com/math/binary-exponentiation/

代码模板：

ll qp(ll a,ll b,ll m){
    ll res = 1;
    a %= m;
    while(b > 0){
        if(b & 1)res = res * a % m;
        a = a * a % m;
        b >>= 1;
    }
    return res;
}

解题思路

就是算\(A^n\)，方法很简单，把快速幂算法中的乘法改成矩阵的乘法就可以了。

AC代码

#include <bits/stdc++.h>
#define debug(a) cout << #a << "=" << a << '\n';
#define il inline
#define int long long

using namespace std;
using ll = long long;
using ull = unsigned long long;

const int N = 110;

int a[N][N], res[N][N], tmp[N][N], m, n,b;

void mx(int a[N][N], int b[N][N], int n) {
    memset(tmp, 0, sizeof(tmp));
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            for (int k = 1; k <= n; k++) {
                tmp[i][j] = (tmp[i][j] + (1LL * a[i][k] * b[k][j]) % m) % m; // 使用 1LL 以防止溢出
            }
        }
    }
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            a[i][j] = tmp[i][j];
        }
    }
}

void qp(int a[N][N], int b, int m) {
    memset(res, 0, sizeof(res));
    for (int i = 1; i <= n; i++) res[i][i] = 1;
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            a[i][j] %= m;
        }
    }
    while (b > 0) {
        if (b & 1) mx(res, a, n);
        mx(a, a, n);
        b >>= 1;
    }
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            a[i][j] = res[i][j];
        }
    }
}

signed main() {
    ios::sync_with_stdio(0), cout.tie(0), cin.tie(0);
    cin >> n >> b;
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            cin >> a[i][j];
        }
    }
    m = 1e9 + 7;
    qp(a, b, m);
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            cout << a[i][j] << ' ';
        }
        cout << '\n';
    }
    return 0;
}

矩阵快速幂加速递推的思路

所谓矩阵加速递推实际上是使用矩阵快速幂这一快速计算矩阵幂次的工具来简化运算，并且简化的是求有彼此有相关关系的数列的运算，如斐波那契额数列的第n项求解（且n往往极大，按照数列表面的数学关系直接递推计算几乎不可能得到）等等；计算这样的递推次数很多的式子，我们需要一个更简化更快捷的计算方式，也就是矩阵快速幂。

首先我们需要想办法把实际问题转换为矩阵的形式来表示，这也是矩阵加速递推最关键的步骤。

例题：斐波那契数列。

其递推关系是：

\(F_1=F_2=1;F_i=F_{i-1}+F_{i-2}(i\geq 3)\)

根据其递推关系我们考虑将矩阵初始值设置为一个行向量（也可以理解成1x2的矩阵，其实列向量还是行向量都可以，只要转移矩阵也是对应的即可）。

\((F_i \qquad F_{i-1})\)

构造一个\(2 \times 2\)的转移矩阵。

我们想要通过乘以转移矩阵，得到数列的下一项，也就是说：

\(\begin{array}{c} (F_i\qquad F_{i-1}) \times \begin{pmatrix} a\qquad b\\ c\qquad d \end{pmatrix} =(F_{i+1}\qquad F_i) \end{array}\)

我们把矩阵乘法展开，得到，

\(\begin{array}{c} a\cdot F_i +c\cdot F_{i-1}=F_{i+1}\\ b\cdot F_i +d\cdot F_{i-1}=F_{i} \end{array}\)

根据数列原本的关系\(F_{i+1}=F_i+F_{i-1}\)，我们可以显而易见求得转移矩阵\(\begin{pmatrix} 1\qquad 1\\1\qquad 0 \end{pmatrix}\)

那么根据转移矩阵表现出来的与原数列的关系，我们可以这样来看，

\(\begin{array}{c} (F_n\qquad F_{n-1}) = (F_2\qquad F_1)\times \begin{pmatrix} 1\qquad 1\\1\qquad 0 \end{pmatrix}^{n-2} \end{array}\)

由此我们就可以解决这个问题了。