矩阵快速幂

矩阵乘法

限制条件 ：\(A\) 的列数等于 \(B\) 的行数

方法：

举个栗子：

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

我们还可以发现，若 \(A\) 是 \(n \times r\) 的矩阵， \(B\) 是 \(r \times m\) 的矩阵，那么 \(C\) 就是 \(n \times m\) 的矩阵

下面给出一个行数、列数都相等的矩阵乘法代码：

inline void Mul(int x[N][N], int y[N][N]) {
    for (int i = 1; i <= n; ++i)
        for (int j = 1; j <= n; ++j)
            ans[i][j] = 0;

    for (int i = 1; i <= n; ++i)
        for (int j = 1; j <= n; ++j)
            for (int k = 1; k <= n; ++k)
                ans[i][j] += x[i][k] * y[k][j];
}

单位矩阵

主对角线上的元素都为 \(1\) ，其余元素都为 \(0\) 的 \(n\) 阶矩阵称作 \(n\) 阶单位矩阵

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

单位矩阵的一个性质：单位矩阵的任意次方都等于它本身。

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

矩阵快速幂

限制条件：矩阵是一个长和宽相等的矩阵

定义：将方阵 \(A\) 自乘 \(n\) 次，记作 \(A^n\)。

因为矩阵乘法满足结合律，那么我们就可以利用快速幂的思想优化

P3390 【模板】矩阵快速幂

如代码所示：

#include <cstdio>
typedef long long ll;
using namespace std;
const ll Mod = 1e9 + 7;
const int N = 1e2 + 7;

template<int SIZE>
struct Matrix {
    ll a[SIZE][SIZE];
    int Row, Col;
    inline Matrix(int r = 0, int c = 0) {
        Row = r, Col = c;

        for (int i = 1; i <= Row; ++i)
            for (int j = 1; j <= Col; ++j)
                a[i][j] = 0;
    }
};

ll k;
int n;

inline Matrix<N> Mul(Matrix<N> &x, Matrix<N> &y) {
    Matrix<N> z(x.Row, y.Col);

    for (int k = 1; k <= x.Col; ++k)
        for (int i = 1; i <= x.Row; ++i)
            for (int j = 1; j <= y.Col; ++j)
                z.a[i][j] = (z.a[i][j] + x.a[i][k] * y.a[k][j] % Mod) % Mod;

    return z;
}

inline Matrix<N> Pow(Matrix<N> p, ll k) {
    Matrix<N> res = p, tmp = p;

    for (--k; k; tmp = Mul(tmp, tmp), k >>= 1)
        if (k & 1)
            res = Mul(res, tmp);

    return res;
}
signed main() {
    scanf("%d%lld", &n, &k);
    Matrix<N> ans(n, n);

    for (int i = 1; i <= n; ++i)
        for (int j = 1; j <= n; ++j)
            scanf("%lld", &ans.a[i][j]);

    ans = Pow(ans, k);

    for (int i = 1; i <= n; ++i) {
        for (int j = 1; j <= n; ++j)
            printf("%lld ", ans.a[i][j]);

        putchar('\n');
    }

    return 0;
}

P1962 斐波那契数列

我们设答案矩阵为 $ F_n = \begin{bmatrix} f_n & f_{n-1} \end{bmatrix}$ ，构建一个初始矩阵 \(A\) 与转移矩阵 \(B\) ，我们希望 \(F_n = F_{n-1} \times B = \begin{bmatrix} f_{n-1} & f_{n-2} \end{bmatrix} \times B\) 推出

因为 \(f_n = f_{n-1} \times 1 + f_{n-2} \times 1\) ，所以 \(B\) 的第一列为 \(\begin{bmatrix} 1 \\ 1 \end{bmatrix}\)

因为 \(f_{n-1} = f_{n-1} \times 1 + f_{n-2} \times 0\) ，所以 \(B\) 的第二列为 \(\begin{bmatrix} 1 \\ 0\end{bmatrix}\)

因此转移矩阵 \(B = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}\)

因为 \(f_1 = f_2 =1\) ，所以我们设初始矩阵 \(A = F_2 = \begin{bmatrix} 1 & 1 \end{bmatrix}\)

综上所述，我们可以推出 \(F_n = A \times B^{n-2}\)

Code：

#include <cstdio>
typedef long long ll;
using namespace std;
const ll Mod = 1e9 + 7;
const int N = 7;

template<int SIZE>
struct Matrix {
    ll a[SIZE][SIZE];
    int Row, Col;
    inline Matrix(int r = 0, int c = 0) {
        Row = r, Col = c;

        for (int i = 1; i <= Row; ++i)
            for (int j = 1; j <= Col; ++j)
                a[i][j] = 0;
    }
};

ll n;

inline Matrix<N> Mul(Matrix<N> &x, Matrix<N> &y) {
    Matrix<N> z(x.Row, y.Col);

    for (int k = 1; k <= x.Col; ++k)
        for (int i = 1; i <= x.Row; ++i)
            for (int j = 1; j <= y.Col; ++j)
                z.a[i][j] = (z.a[i][j] + x.a[i][k] * y.a[k][j] % Mod) % Mod;

    return z;
}

inline Matrix<N> Pow(Matrix<N> p, ll k) {
    Matrix<N> res = p, tmp = p;

    for (--k; k; tmp = Mul(tmp, tmp), k >>= 1)
        if (k & 1)
            res = Mul(res, tmp);

    return res;
}
signed main() {
    scanf("%lld", &n);

    if (n <= 2) {
        putchar('1');
        return 0;
    }

    Matrix<N> A(1, 2), B(2, 2), ans(1, 2);
    A.a[1][1] = A.a[1][2] = 1;
    B.a[1][1] = B.a[1][2] = B.a[2][1] = 1, B.a[2][2] = 0;
    B = Pow(B, n - 2), ans = Mul(A, B);
    printf("%lld", ans.a[1][1]);
    return 0;
}

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

我们设答案矩阵为 $ F_n = \begin{bmatrix} f_n & f_{n-1} & f_{n-2} \end{bmatrix}$ ，构建一个初始矩阵 \(A\) 与转移矩阵 \(B\) ，我们希望 \(F_n = F_{n-1} \times B = \begin{bmatrix} f_{n-1} & f_{n-2} & f_{n-3} \end{bmatrix} \times B\) 推出

因为 \(f_n = f_{n-1} \times 1 + f_{n-2} \times 0 + f_{n-3} \times 1\) ，所以 \(B\) 的第一列为 \(\begin{bmatrix} 1 \\ 0 \\ 1 \end{bmatrix}\)

因为 $f_{n-1} = f_{n-1} \times 1 + f_{n-2} \times 0 + f_{n-3} \times 0 $ ，所以 \(B\) 的第二列为 \(\begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}\)

因为 $f_{n-2} = f_{n-1} \times 0 + f_{n-2} \times 1 + f_{n-3} \times 0 $ ，所以 \(B\) 的第二列为 \(\begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix}\)

因此转移矩阵 \(B = \begin{bmatrix} 1 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{bmatrix}\)

因为 \(f_1 = f_2 = f_3 = 1\) ，所以我们设初始矩阵 \(A = F_3 = \begin{bmatrix} 1 & 1 & 1 \end{bmatrix}\)

综上所述，我们可以推出 \(F_n = A \times B^{n-3}\)

Code：

#include <cstdio>
typedef long long ll;
using namespace std;
const ll Mod = 1e9 + 7;
const int N = 7;

template<int SIZE>
struct Matrix {
    ll a[SIZE][SIZE];
    int Row, Col;
    inline Matrix(int r = 0, int c = 0) {
        Row = r, Col = c;

        for (int i = 1; i <= Row; ++i)
            for (int j = 1; j <= Col; ++j)
                a[i][j] = 0;
    }
};

ll T, n;

inline Matrix<N> Mul(Matrix<N> &x, Matrix<N> &y) {
    Matrix<N> z(x.Row, y.Col);

    for (int k = 1; k <= x.Col; ++k)
        for (int i = 1; i <= x.Row; ++i)
            for (int j = 1; j <= y.Col; ++j)
                z.a[i][j] = (z.a[i][j] + x.a[i][k] * y.a[k][j] % Mod) % Mod;

    return z;
}

inline Matrix<N> Pow(Matrix<N> p, ll k) {
    Matrix<N> res = p, tmp = p;

    for (--k; k; tmp = Mul(tmp, tmp), k >>= 1)
        if (k & 1)
            res = Mul(res, tmp);

    return res;
}
signed main() {
    scanf("%lld", &T);

    for (; T; --T) {
        scanf("%lld", &n);

        if (n <= 3) {
            puts("1");
            continue;
        }

        Matrix<N> A(1, 3), B(3, 3), ans(1, 2);
        A.a[1][1] = 1, A.a[1][2] = 1, A.a[1][3] = 1;
        B.a[1][1] = 1, B.a[1][2] = 1, B.a[1][3] = 0;
        B.a[2][1] = 0, B.a[2][2] = 0, B.a[2][3] = 1;
        B.a[3][1] = 1, B.a[3][2] = 0, B.a[3][3] = 0;
        B = Pow(B, n - 3), ans = Mul(A, B);
        printf("%lld\n", ans.a[1][1]);
    }

    return 0;
}

本栏目推荐文章
• 矩阵乘法代码
• 螺旋矩阵 59
• 数学基础(一)-标量、向量、矩阵、张量以及各范数的含义
• 矩阵相关知识图谱
• 实现自定义算子：逆矩阵inverse
• 机器学习-矩阵
• 顺时针旋转矩阵
• OpenCL使用local内存优化矩阵乘法例子
• 营业日志 —— 矩阵加速
• CF1864H Asterism Stream【概率 DP，矩阵优化】