矩阵快速幂模板

最新推荐文章于 2024-10-24 19:05:02 发布

hesorchen

最新推荐文章于 2024-10-24 19:05:02 发布

阅读量1.2k

点赞数 5

CC 4.0 BY-SA版权

分类专栏： # 矩阵快速幂文章标签：线性代数

本文链接：https://blog.youkuaiyun.com/hesorchen/article/details/119918350

矩阵快速幂专栏收录该内容

1 篇文章

订阅专栏

本文详细介绍了矩阵乘法的基本概念，包括矩阵乘法的定义、性质，以及单位矩阵和零矩阵。接着讲解了矩阵快速幂的原理，通过C++代码展示了如何实现矩阵乘法和矩阵快速幂，并提供了将递推式转化为矩阵乘法的模板。最后，通过实例展示了如何使用矩阵快速幂解决斐波那契数列问题，进一步说明了其在算法加速中的应用价值。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

前置知识

矩阵

矩阵是一个 $\times m$ 的阵列，下面是一个 $3\times 3$ 的矩阵：

$\begin{bmatrix} 1&0&0\\\\ 0&1&0\\\\ 0&0&1\\ \end{bmatrix}$

矩阵乘法

若矩阵A的大小为 $n\times m$ ，另一矩阵B的大小为 $m\times p$ ，则两个矩阵可以做乘法，得到的矩阵C大小为 $n * p$ 。具体的乘法规则如下：

矩阵A：
$\begin{bmatrix} a_{11}&a_{12}&a_{13}\\\\ a_{21}&a_{22}&a_{23}\\ \end{bmatrix}$

矩阵B：
$\begin{bmatrix} b_{11}&b_{12}\\\\ b_{21}&b_{22}\\\\ b_{31}&b_{32}\\ \end{bmatrix}$

相乘得到矩阵C：
$\begin{bmatrix} a_{11}\times b_{11}+a_{12}\times b_{21}+a_{13}\times b_{31} &a_{11}\times b_{12}+a_{12}\times b_{22}+a_{13}\times b_{32}\\\\ a_{21}\times b_{11}+a_{22}\times b_{21}+a_{23}\times b_{31}&a_{21}\times b_{12}+a_{22}\times b_{22}+a_{23}\times b_{32} \end{bmatrix}$

C++代码实现矩阵乘法:

for (int i = 1; i <= n; i++)
    for (int j = 1; j <= n; j++)
        for (int k = 1; k <= n; k++)
            C[i][j] += A[i][k] * B[k][j];

性质：

矩阵乘法不满足交换律，但是满足结合律。即 $(A B) C = A (B C)$ 。
单位矩阵：单位矩阵起着特殊的作用，如同数的乘法中的1，这种矩阵被称为单位矩阵。它是个方阵，从左上角到右下角的对角线（称为主对角线）上的元素均为1。除此以外全都为0。根据单位矩阵的特点，任何矩阵与单位矩阵相乘都等于本身。
零矩阵：零矩阵即所有元素皆为0的矩阵。

矩阵快速幂

类似于快速幂用于加速乘法，矩阵快速幂用于加速矩阵乘法。我们可以将 $O (n)$ 的乘法优化到 $O (l o g n)$ 。

可以对比一下快速幂与矩阵快速幂：

long long quick_pow(long long a, long long b)//快速幂
{
    long long res = 1; //1
    while (b)
    {
        if (b & 1)
            res = res * a;
        b /= 2;
        a = a * a;
    }
    return res;
}

Matrix quick_multi(Matrix a, int t) //矩阵快速幂
{
    Matrix res; //单位矩阵
    for (int i = 1; i <= N; i++)
        for (int j = 1; j <= N; j++)
            res.M[i][j] = (i == j);
    while (t)
    {
        if (t & 1)
            res = multi(res, a);
        t /= 2;
        a = multi(a, a);
    }
    return res;
}

应用

如果我们能够将一个递推式表示成矩阵乘法的形式，那么就可以用矩阵快速幂加速。例如经典的矩阵快速幂加速斐波那契数列。

我们记斐波那契第 $i$ 项为 $F_i$ ，可以将 $F_i=F_{i-1}+F_{i-2}$ 写成矩阵乘法的形式：

$\begin{bmatrix} F_{n}&F_{n-1}\\ \end{bmatrix}=\begin{bmatrix} F_{n-1}&F_{n-2}\\ \end{bmatrix}\times \begin{bmatrix} 1&1\\ 1&0\\ \end{bmatrix}$

更一般的，可以写作

$\begin{bmatrix} F_{n}&F_{n-1}\\ \end{bmatrix}=\begin{bmatrix} F_{2}&F_{1}\\ \end{bmatrix}\times \begin{bmatrix} 1&1\\ 1&0\\ \end{bmatrix}^{n-2}$

这样就可以用 $O (l o g n)$ 的时间复杂度得到 $F_n$

套路

如果已经得到了矩阵乘法，那矩阵快速幂就是一个模板，一般的矩阵快速幂题目难点在于写出递推式以及将递推式转化成矩阵乘法。递推式因题而异，而将递推式转换为矩阵乘法有着一个比较常用的套路，例如要将 $F_n=F_{n-1}+n^3+n^2+1$ 转化为矩阵乘法。
我们先写出 $F_n$ 以及 $F_{n-1}$ 的表达式：

$F_n=F_{n-1}+n^3+n^2+1$
$F_{n-1}=F_{n-2}+(n-1)^3+(n-1)^2+1$

第一列一般就是 $F_n$ 的递推式，最后一列是 $F_{n-1}$ 的递推式，这样才能从 $F_{n-1}$ 递推到 $F_{n}$ ：

$\begin{bmatrix} F_{n}\\\\ F_{n-1}\\\\ n^3\\\\ n^2\\\\ 1\\ \end{bmatrix} = \begin{bmatrix} & & & & \\\\ & & & & \\\\ & & & & \\\\ & & & & \\\\ & & & & & &&& \\ \end{bmatrix} \times \begin{bmatrix} F_{n-1}\\\\ F_{n-2}\\\\ (n-1)^3\\\\ (n-1)^2\\\\ 1\\ \end{bmatrix}$

其他项的原则是缺什么补什么，当要表示 $n^3$ 以及 $n^2$ 时，我们会发现，需要补一个 $n$ 才能表示出它们：

$\begin{bmatrix} F_{n}\\\\ F_{n-1}\\\\ n^3\\\\ n^2\\\\ n\\\\ 1\\ \end{bmatrix} = \begin{bmatrix} & & & & \\\\ & & & & \\\\ & & & & \\\\ & & & & \\\\ & & & & \\\\ & & & & & &&& \\ \end{bmatrix} \times \begin{bmatrix} F_{n-1}\\\\ F_{n-2}\\\\ (n-1)^3\\\\ (n-1)^2\\\\ n\\\\ 1\\ \end{bmatrix}$

需要的项都够了之后，只需仔细的写出系数矩阵即可。

$\begin{bmatrix} F_{n}\\\\ F_{n-1}\\\\ n^3\\\\ n^2\\\\ n\\\\ 1\\ \end{bmatrix} = \begin{bmatrix} 1&0&1&4&5&-2\\\\ 1&0&0&0&0&0\\\\ 0&0 &1 &3 &3 &-2\\\\ 0& 0& 0& 1& 2&-1\\\\ 0& 0& 0& 0& 1&0\\\\ 0& 0& 0& 0& 0&1\\ \end{bmatrix} \times \begin{bmatrix} F_{n-1}\\\\ F_{n-2}\\\\ (n-1)^3\\\\ (n-1)^2\\\\ n-1\\\\ 1\\ \end{bmatrix}$

代码

/*
Fn = Fn-1 + n^3 + n^2 + 1
*/
#include <bits/stdc++.h>
using namespace std;
const int N = 8;
const int mod = 13331;
int LEN;
struct Matrix
{
    int M[N][N];
};
void init(Matrix &a, int opt) //opt=0 初始化为零矩阵
{
    for (int i = 1; i <= LEN; i++)
        for (int j = 1; j <= LEN; j++)
            a.M[i][j] = (opt == 1 && i == j);
}
Matrix multi(Matrix a, Matrix b)
{
    Matrix res;
    init(res, 0);
    for (int i = 1; i <= LEN; i++)
        for (int j = 1; j <= LEN; j++)
            for (int k = 1; k <= LEN; k++)
            {
                res.M[i][j] += a.M[i][k] * b.M[k][j];
                res.M[i][j] %= mod;
            }
    return res;
}
Matrix quick_multi(Matrix a, int t) //矩阵快速幂
{
    Matrix res;
    init(res, 1); //初始化为单位矩阵
    while (t)
    {
        if (t & 1)
            res = multi(res, a);
        t /= 2;
        a = multi(a, a);
    }
    return res;
}
int BEG[7][7] = {{},
                 {0, 1},
                 {0, 0},
                 {0, 1},
                 {0, 1},
                 {0, 1},
                 {0, 1}};
int BASE[7][7] = {
    {},
    {0, 1, 0, 1, 4, 5, 3},
    {0, 1, 0, 0, 0, 0, 0},
    {0, 0, 0, 1, 3, 3, 1},
    {0, 0, 0, 0, 1, 2, 1},
    {0, 0, 0, 0, 0, 1, 1},
    {0, 0, 0, 0, 0, 0, 1}};
void solve()
{
    int n;
    cin >> n;
    LEN = 6;
    Matrix ans, beg, base;
    for (int i = 1; i <= LEN; i++)
        for (int j = 1; j <= LEN; j++)
            base.M[i][j] = BASE[i][j], beg.M[i][j] = BEG[i][j];
    base = quick_multi(base, n - 1);
    init(ans, 0);
    for (int i = 1; i <= LEN; i++)
        for (int j = 1; j <= 1; j++)
            for (int k = 1; k <= LEN; k++)
            {
                ans.M[i][j] += base.M[i][k] * beg.M[k][j];
                ans.M[i][j] %= mod;
            }
    cout << ans.M[1][1] << endl;
}
int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    int t;
    cin >> t;
    while (t--)
        solve();
    return 0;
}

例题

POJ3070 Fibonacci

代码

#include <iostream>
using namespace std;

const int N = 1e5 + 5;
const int mod = 10000;
struct Matrix
{
    int M[5][5];
};

Matrix multi(Matrix a, Matrix b)
{
    Matrix ans;
    for (int i = 1; i <= 2; i++)
        for (int j = 1; j <= 2; j++)
            ans.M[i][j] = 0;
    for (int i = 1; i <= 2; i++)
        for (int j = 1; j <= 2; j++)
            for (int k = 1; k <= 2; k++)
            {
                ans.M[i][j] += a.M[i][k] * b.M[k][j];
                ans.M[i][j] %= mod;
            }
    return ans;
}

Matrix quick_multi(Matrix base, int t)
{
    Matrix res;
    res.M[1][1] = res.M[2][2] = 1;
    res.M[1][2] = res.M[2][1] = 0;
    while (t)
    {
        if (t & 1)
            res = multi(res, base);
        t /= 2;
        base = multi(base, base);
    }
    return res;
}
void solve()
{
    int n;
    while (cin >> n)
    {
        if (n <= 2)
        {
            if (n == -1)
                return;
            if (n == 0)
                cout << 0 << endl;
            else
                cout << 1 << endl;
            continue;
        }
        Matrix beg;
        beg.M[1][1] = 1;
        beg.M[1][2] = 1;
        Matrix base;
        base.M[1][1] = 0;
        base.M[1][2] = 1;
        base.M[2][1] = 1;
        base.M[2][2] = 1;
        base = quick_multi(base, n - 2);
        Matrix ans;
        for (int i = 1; i <= 2; i++)
            for (int j = 1; j <= 2; j++)
                ans.M[i][j] = 0;
        for (int i = 1; i <= 1; i++)
            for (int j = 1; j <= 2; j++)
                for (int k = 1; k <= 2; k++)
                {
                    ans.M[i][j] += beg.M[i][k] * base.M[k][j];
                    ans.M[i][j] %= mod;
                }
        cout << ans.M[1][2] << endl;
    }
}
int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0);
    solve();
    return 0;
}
/*


*/