【算法分享】快速幂与矩阵快速幂算法详解

在程序竞赛和实际开发中，我们经常遇到需要计算\(a^n\)这样的幂运算。当n很大时（如\(10^9\)），直接用循环运算肯定会TLE。这时候就需要用到快速幂算法，而当我们需要计算矩阵的幂时，矩阵快速幂就派上用场了。

🔥 第一部分：普通快速幂

🤔 问题引入

假设我们要计算 \(2^{10}\)，你会怎么算？

方法1：暴力计算 😫

ll b=2;
for(int i=1;i<10;i++)
{
  b*=b;
}

2^10 = 2 × 2 × 2 × 2 × 2 × 2 × 2 × 2 × 2 × 2

需要进行9次乘法运算。

方法2：聪明的计算 🧠✨

2^1 = 2
2^2 = 4
2^4 = 16
2^8 = 256
2^10 = 2^8 × 2^2 = 256 × 4 = 1024

只需要4次乘法运算！🎉

二进制分解的奥秘

关键在于理解：任何数都可以用二进制表示！

比如 \(10=(1010)_{2}=1×2^3+0×2^2+1×2^1+0×2^0=8+2\)

所以：\(a^{10}=a^{(8+2)}=a^8×a^2\)

一般化：如果$ n={(b_{k} b_{k-1}...b_{1}b_{0})}_{2}$，那么：

\(a^n = a^{(b_k×2^k + b_(k-1)×2^(k-1) + ... + b_1×2^1 + b_0×2^0)} \)

\(= \prod\limits_{i=0}^ka^{(b_i×2^i)}\)

其中\(b_{i}​∈{0,1}，当b_{i}=1\) 时才参与乘积。

📋 算法步骤详解

让我们用 \(3^{13}\) 来演示整个过程：

1：将指数转为二进制

\(13=(1101)_{2}\)

2：从右到左处理每一位 👆

初始: result = 1, base = 3, n = 13 = (1101)₂

第1次循环: n的最后一位是1 
- 因为最后一位是1，所以 result = result × base = 1 × 3 = 3
- base = base × base = 3 × 3 = 9
- n = n >> 1 = (110)₂ = 6

第2次循环: n的最后一位是0 
- 因为最后一位是0，所以 result 不变 = 3
- base = base × base = 9 × 9 = 81
- n = n >> 1 = (11)₂ = 3

第3次循环: n的最后一位是1 ✅
- 因为最后一位是1，所以 result = result × base = 3 × 81 = 243
- base = base × base = 81 × 81 = 6561  
- n = n >> 1 = (1)₂ = 1

第4次循环: n的最后一位是1 ✅
- 因为最后一位是1，所以 result = result × base = 243 × 6561 = 1594323
- base = base × base = 6561 × 6561 (不需要了)
- n = n >> 1 = 0

循环结束，返回 result = 1594323 🎯

代码实现

ll ksm1(ll base ,ll n) {
    ll res=1;
    while(n>0) {
        if(n&1) {
            res*=base;
        }
        base *= base;
        n>>=1;
    }
    return res;
}


//带模运算(防止溢出)
ll ksm2(ll base ,ll n,ll m) {
    ll res=1;
    while (n>0) {
        if (n&1) {
            res=res*base%m;
        }
        base=base*base%m;
        n>>=1;
    }
    return res;
}

1.5 📊 复杂度分析

时间复杂度：O(logn)，因为每次循环指数都除以2 ⚡

空间复杂度：O(1)，只使用了常数个变量 💾

优势：将原本 O(n) 的暴力算法优化到 O(logn) 🚀

📐 第二部分：矩阵快速幂

在学习矩阵快速幂之前，我们了解一些矩阵的基本概念。

矩阵乘法

重要规则：只有当第一个矩阵的列数等于第二个矩阵的行数时，才能相乘。⚠️

设 A 是$ m×n $矩阵，B 是 $ n×p$ 矩阵，则 C\(=A×B\) 是 \(m×p\) 矩阵，其中： \(C[i][j]\sum\limits_{k=1}^n​A[i][k]×B[k][j] \)

例子

\(\begin{bmatrix} 1 & 2\\ 3 & 4 \end{bmatrix}\cdot\begin{bmatrix}5 & 6\\7 & 8\end{bmatrix}=\begin{bmatrix}1\cdot5+2\cdot7 & 1\cdot6+2\cdot8\\3\cdot5+4\cdot7 & 3\cdot6+4\cdot8\end{bmatrix}=\begin{bmatrix}19 & 22\\43 & 50\end{bmatrix} \)

单位矩阵

单位矩阵是矩阵中的"1"，对角线为1，其他位置为0​​

\(I_2 = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}\)

\(I_3 = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}\)

假设我们要计算一个矩阵的高次幂，比如：

\(A^{1000000} =\underbrace{A \times A \times \cdots \times A}_{\text{1000000 次}}\)

如果用暴力方法，需要999999次矩阵乘法，每次乘法是$ O(n^3)$，总复杂度是 \(O(n^3×10^6)\)，显然会超时！😱

这时候我们就需要矩阵快速幂！🦸‍♂️

💡 核心思想

矩阵快速幂的思想与普通快速幂完全相同：

将指数进行二进制分解 🔢

利用矩阵乘法的结合律 🔗

通过平方来快速计算 ⚡

📝 算法步骤

以计算 \(A^5\) 为例，\(5=(101)_{2}\)：

初始: result = I（单位矩阵）, base = A, exp = 5

第1次循环: exp的最后一位是1 
- result = result × base = I × A = A  
- base = base × base = A × A = A²
- exp = exp >> 1 = 2
  
第2次循环: exp的最后一位是0 
- result 不变 = A
- base = base × base = A² × A² = A⁴  
- exp = exp >> 1 = 1

第3次循环: exp的最后一位是1 
- result = result × base = A × A⁴ = A⁵
- exp = exp >> 1 = 0

循环结束，返回 A⁵ 🎉

代码实现

// 定义一个矩阵结构体 Me
struct Me {
    ll a[105][105];  // 存放矩阵元素，下标范围 1..n
    
    // 矩阵乘法运算符重载：C = A * B
    Me operator*(const Me &b) {
        Me res;  // 结果矩阵
        for (int i = 1; i <= n; i++) {       // 枚举结果矩阵的行
            for (int j = 1; j <= n; j++) {   // 枚举结果矩阵的列
                res.a[i][j] = 0;             // 初始化结果矩阵的元素为 0
                for (int k = 1; k <= n; k++) {   
                    // res[i][j] = Σ (a[i][k] * b[k][j])
                    // 这里取模 mm，防止溢出
                    res.a[i][j] = (a[i][k] * b.a[k][j] % mm + res.a[i][j] % mm) % mm;
                }
            }
        }
        return res;  // 返回结果矩阵
    }
} a; // 定义一个全局矩阵 a


Me ksm(Me k, ll x) {
    Me res;
    // 初始化单位矩阵
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            res.a[i][j] = (i == j) ? 1 : 0;
        }
    }
    while (x) {
        if (x&1)
            res=res*k;
        k=k*k;
        x>>=1;
    }
    return res;
}

3.5 📊 复杂度分析

时间复杂度：\(O(\)\(n^3logk\))，其中n是矩阵大小，k是指数 ⏱️

空间复杂度：\(O(n^2)\) 💾

优势：将$ O$$(n^3k)$ 优化到$ O(n^3logk)$

实战题目

5.1 📝 洛谷P3390 【模板】矩阵快速幂

题目：给定n×n矩阵A，求Ak。

分析：这是矩阵快速幂的模板题，直接套用算法即可。✅

#include <iostream>
#include <queue>
#include <cstring>
using namespace std;
typedef unsigned long long ll;

ll n,m,k;
const ll mm=1e9+7;
struct Me{
    ll a[105][105];
    Me operator*(const Me &b) {
        Me res;
        for (int i=1;i<=n;i++) {
            for (int j=1;j<=n;j++) {
                res.a[i][j]=0;
                for (int k=1;k<=n;k++) {
                    res.a[i][j]=(a[i][k]*b.a[k][j]%mm+res.a[i][j]%mm)%mm;
                }
            }
        }
        return res;
    }
}a;

Me ksm(Me k, ll x) {
    Me res;
    // 初始化单位矩阵
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            res.a[i][j] = (i == j) ? 1 : 0;
        }
    }
    while (x) {
        if (x&1)
            res=res*k;
        k=k*k;
        x>>=1;
    }
    return res;
}


int main() {
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);

    ll b=2;
    for(int i=1;i<10;i++)
    {
        b*=b;
        cout<<b<<endl;
    }
    cout<<b<<endl;
    return 0;
}