矩阵快速幂模板

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

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

前置知识

矩阵

矩阵是一个 n × m n \times m n×m的阵列,下面是一个 3 × 3 3\times 3 3×3的矩阵:

[ 1 0 0 0 1 0 0 0 1 ] \begin{bmatrix} 1&0&0\\\\ 0&1&0\\\\ 0&0&1\\ \end{bmatrix} 100010001

矩阵乘法

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

矩阵A:
[ a 11 a 12 a 13 a 21 a 22 a 23 ] \begin{bmatrix} a_{11}&a_{12}&a_{13}\\\\ a_{21}&a_{22}&a_{23}\\ \end{bmatrix} a11a21a12a22a13a23

矩阵B:
[ b 11 b 12 b 21 b 22 b 31 b 32 ] \begin{bmatrix} b_{11}&b_{12}\\\\ b_{21}&b_{22}\\\\ b_{31}&b_{32}\\ \end{bmatrix} b11b21b31b12b22b32

相乘得到矩阵C:
[ a 11 × b 11 + a 12 × b 21 + a 13 × b 31 a 11 × b 12 + a 12 × b 22 + a 13 × b 32 a 21 × b 11 + a 22 × b 21 + a 23 × b 31 a 21 × b 12 + a 22 × b 22 + a 23 × b 32 ] \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} a11×b11+a12×b21+a13×b31a21×b11+a22×b21+a23×b31a11×b12+a12×b22+a13×b32a21×b12+a22×b22+a23×b32

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];

性质

  1. 矩阵乘法不满足交换律,但是满足结合律。即 ( A B ) C = A ( B C ) (AB)C=A(BC) (AB)C=A(BC)

  2. 单位矩阵:单位矩阵起着特殊的作用,如同数的乘法中的1,这种矩阵被称为单位矩阵。它是个方阵,从左上角到右下角的对角线(称为主对角线)上的元素均为1。除此以外全都为0。根据单位矩阵的特点,任何矩阵与单位矩阵相乘都等于本身。

  3. 零矩阵:零矩阵即所有元素皆为0的矩阵。

矩阵快速幂

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

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

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 i i项为 F i F_i Fi,可以将 F i = F i − 1 + F i − 2 F_i=F_{i-1}+F_{i-2} Fi=Fi1+Fi2写成矩阵乘法的形式:

[ F n F n − 1 ] = [ F n − 1 F n − 2 ] × [ 1 1 1 0 ] \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} [FnFn1]=[Fn1Fn2]×[1110]

更一般的,可以写作

[ F n F n − 1 ] = [ F 2 F 1 ] × [ 1 1 1 0 ] n − 2 \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} [FnFn1]=[F2F1]×[1110]n2

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

套路

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

F n = F n − 1 + n 3 + n 2 + 1 F_n=F_{n-1}+n^3+n^2+1 Fn=Fn1+n3+n2+1
F n − 1 = F n − 2 + ( n − 1 ) 3 + ( n − 1 ) 2 + 1 F_{n-1}=F_{n-2}+(n-1)^3+(n-1)^2+1 Fn1=Fn2+(n1)3+(n1)2+1

第一列一般就是 F n F_n Fn的递推式,最后一列是 F n − 1 F_{n-1} Fn1的递推式,这样才能从 F n − 1 F_{n-1} Fn1递推到 F n F_{n} Fn

[ F n F n − 1 n 3 n 2 1 ] = [ ] × [ F n − 1 F n − 2 ( n − 1 ) 3 ( n − 1 ) 2 1 ] \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} FnFn1n3n21=×Fn1Fn2(n1)3(n1)21

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

[ F n F n − 1 n 3 n 2 n 1 ] = [ ] × [ F n − 1 F n − 2 ( n − 1 ) 3 ( n − 1 ) 2 n 1 ] \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} FnFn1n3n2n1=×Fn1Fn2(n1)3(n1)2n1

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

[ F n F n − 1 n 3 n 2 n 1 ] = [ 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 ] × [ F n − 1 F n − 2 ( n − 1 ) 3 ( n − 1 ) 2 n − 1 1 ] \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} FnFn1n3n2n1=110000000000101000403100503210202101×Fn1Fn2(n1)3(n1)2n11

代码

/*
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;
}
/*


*/

参考资料

  1. 矩阵快速幂模板

  2. 矩阵快速幂(原理+模板)

### 矩阵快速幂算法的实现 矩阵快速幂是一种高效的算法,用于计算矩阵的高次幂。它基于分治的思想以及矩阵乘法的结合律来降低时间复杂度。以下是矩阵快速幂的一个通用代码模板: #### Python 实现 ```python import numpy as np def matrix_multiply(A, B, mod=None): """矩阵相乘""" rows_A, cols_A = len(A), len(A[0]) rows_B, cols_B = len(B), len(B[0]) if cols_A != rows_B: raise ValueError("无法进行矩阵乘法") result = [[0 for _ in range(cols_B)] for __ in range(rows_A)] for i in range(rows_A): for j in range(cols_B): temp_sum = 0 for k in range(cols_A): temp_sum += A[i][k] * B[k][j] if mod is not None: temp_sum %= mod result[i][j] = temp_sum return result def matrix_power(matrix, n, mod=None): """矩阵快速幂""" size = len(matrix) identity_matrix = [[int(i == j) for j in range(size)] for i in range(size)] result = identity_matrix base = matrix while n > 0: if n % 2 == 1: result = matrix_multiply(result, base, mod=mod) base = matrix_multiply(base, base, mod=mod) n //= 2 return result ``` 上述代码实现了两个核心函数: - `matrix_multiply`:完成两个矩阵之间的乘法操作,并支持模运算[^1]。 - `matrix_power`:通过快速幂的方式高效地计算矩阵的高次幂。 #### C++ 实现 对于更注重性能的语言如C++,也可以提供类似的实现方式: ```cpp #include <vector> using namespace std; // 定义矩阵大小和取模值 const int MOD = 1e9 + 7; typedef vector<vector<long long>> Matrix; Matrix multiply(const Matrix &A, const Matrix &B){ int r = A.size(), c = B[0].size(); Matrix C(r, vector<long long>(c, 0)); for(int i = 0;i < r;i++) { for(int j = 0;j < c;j++) { for(int k = 0;k < (int)B.size();k++) { C[i][j] = (C[i][j] + A[i][k]*B[k][j])%MOD; } } } return C; } Matrix power(Matrix base, long long exp){ int sz = base.size(); Matrix res(sz, vector<long long>(sz, 0)); // 单位矩阵初始化 for(int i = 0;i < sz;i++) res[i][i] = 1; while(exp > 0){ if(exp & 1){ // 如果当前指数为奇数 res = multiply(res, base); } base = multiply(base, base); // 平方更新基底 exp >>= 1; // 右移一位相当于除以2 } return res; } ``` 以上代码同样包含了两部分功能: - `multiply` 函数负责执行矩阵间的乘法并处理大整数溢出问题[^4]。 - `power` 函数则采用快速幂的方法加速矩阵幂次的计算。 #### 应用实例——斐波那契数列 假设我们需要使用矩阵快速幂求解第 \(n\) 项斐波那契数列,则可以通过如下构造矩阵来进行计算: \[ M = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}, V_0 = \begin{bmatrix} F(1)\\ F(0) \end{bmatrix} = \begin{bmatrix} 1\\ 0 \end{bmatrix}. \] 那么有 \(\text{{result}} = M^{n-1} \times V_0\) 表示最终的结果向量[^3]。 --- ###
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

hesorchen

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值