Reference(Aka thanks)(可能"Reference"在这里的用法不太对，欢迎提醒我来纠正)
- Wiki
  - 时间复杂度 - 维基百科，自由的百科全书
  - 快速幂 - OI Wiki
- 基础
  - [[Binary Exponentiation]]
  - [[Linear Algebra|矩阵乘法]]
- 视频
  - 509 Luogu P3390【模板】矩阵快速幂 bilibili
- 题目
  - P1226 【模板】快速幂 - 洛谷 | 计算机科学教育新生态
  - P3390 【模板】矩阵快速幂 - 洛谷 | 计算机科学教育新生态
- 讨论(想法)
  - 似乎还可以通过剪枝优化，下次()
  - 关于快速幂的名字笔者并无好感，更喜欢称其为“二进制求幂”，非常直观

把基础拿出来讲 -- 矩阵乘法

矩阵乘法:

TL;DR: 前一个矩阵逐行 与 后一个矩阵逐列 相乘求和 -- 结果放入结果矩阵 Z 字排列

设 A 为 m*p 的矩阵，B 为 p*n 的矩阵，称 m*n 的矩阵 C 为 A 与 B 的乘积，则 C 的第 i 行第 j 列的元素为

c_{ij} = \sum_{k=1}^p a_{ik}b_{kj} = a_{i1}b_{1j} + a_{i2}b_{2j} + \cdots + a_{ip}b_{pj}

m 为行，p 为列

Complexity
- 朴素求幂复杂度为 \mathcal{O(N)}，根据主定理(<-朋友告诉的，笔者还不会主定理，算法导论，请)，分析得到快速幂的时间复杂度为 \mathcal{O(logN)}，
  实际上还可以通过费马小定理加速
  - if isPrime(m) x^{n \pmod{m - 1}}
- 进而来到矩阵快速幂的复杂度: 设矩阵维度为 m，矩阵乘法的复杂度为 m^3，快速幂的复杂度为 log_n，故矩阵快速幂的时间复杂度为 m^3*log n，而 m 为常数，而根据时间复杂度 - 维基百科，自由的百科全书中所述，在大O表示法中常数不计。矩阵快速幂复杂度即为 \mathcal{O(log N)}

初学编程的伙伴可能对于循环中 i, j 不是很敏感，或者容易弄混，这里特意排出这个矩阵，供参考:

一个 m \times n 的矩阵是一个由 m 行 n 列元素排列成的矩形阵列。即形如

A = \underbrace{ \left. \begin{bmatrix} a_{1 1} & a_{1 2} & \cdots & a_{1 n} \\ a_{2 1} & a_{2 2} & \cdots & a_{2 n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m 1} & a_{m 2} & \cdots & a_{m n} \\ \end{bmatrix} \right\} }_{n} m \text{.}

矩阵快速幂的思路

划到最后看代码，注释写的比较清晰了(自认为)

最开始题目所对应的答案(Template Code)

快速幂

Python

def quick_pow(x, n, mod = None):
    res = 1
    while n:
        if n & 1:
            res = res * x % mod
        x = x * x % mod
        # print(x)
        n >>= 1
    return res

print(f'{quick_pow(2, 4) = }') # 16

C++

#include <iostream>
// 快速幂
long long quick_pow(long long a, long long b, long long p){
    long long ans = 1;
    while (b) {
        if (b & 1) ans = ans * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return ans;
}

int main() {
    long long a, b, p;
    std::cin >> a >> b >> p;
    // 2^10 mod 9=7
    std::cout << a << "^" << b << " mod " << p << "=" << quick_pow(a, b, p) << std::endl;
    return 0;
}

矩阵快速幂

C++
使用 clang 编译此程序需要链接 libc++ 如下:

clang++ -std=c++17 -stdlib=libc++ -o {二进制文件保存为} {待编译的C++源代码}.cpp

#include <iostream> // io
#include <cstring> // memset

const long long MOD = 1e9 + 7; // 10**9 + 7

struct matrix
{
    long long coefficient[105][105];
    matrix() {
        memset(coefficient, 0, sizeof(coefficient)); // 初始化为0
    }
} matrixA, result; // 初始所给的矩阵matrixA 结果矩阵result

int n;
long long k;

// 矩阵乘法 LaTex: 设matrixA为m*p的矩阵 B为p*n的矩阵 则称m*n的矩阵C为matrixA与B的乘积 则C的第i行第j列的元素为
// [imath:0]c_{ij} = \sum_{k=1}^p a_{ik}b_{kj} = a_{i1}b_{1j} + a_{i2}b_{2j} + \cdots + a_{ip}b_{pj}[/imath:0]
// m为行 p为列
// TL;DR: 前一个矩阵逐行 与 后一个矩阵逐列 相乘求和 -- 结果放入结果矩阵Z字排列
matrix operator*(matrix &a, matrix &b) { // 重载乘号 实现矩阵乘法
    matrix ans; // 定义矩阵时会自动调用构造函数 也就是 matrix() 使其初始化为0
    for (int i = 1; i <= n; i++) {
        for (int k = 1; k <= n; k++) {
            for (int j = 1; j <= n; j++) {
                // 实现乘积求和的过程 原来的值 + 乘积 
                ans.coefficient[i][j] = (ans.coefficient[i][j] + a.coefficient[i][k] * b.coefficient[k][j]) % MOD;
                // 第一维下标尽量放在外层循环 第二维下标放在内层循环 这样可以提高效率
                // 也就是 i j k --> i k j
            }
        }
    }
    return ans;
}

void quickPower(long long exponentiation) { // 快速幂
    for (int i = 1; i <= n; i++) { // 构造单位矩阵
        result.coefficient[i][i] = 1; // 主对角线为1
    } // 单位矩阵乘任何矩阵都是原矩阵
    while (exponentiation) {
        if (exponentiation & 1) {
            result = result * matrixA; // c++智能识别变量类型来决定是否调用矩阵乘法函数 使重载通用
        }
        matrixA = matrixA * matrixA;
        exponentiation >>= 1;
    }
}

int main() {
    scanf("%d%lld",&n,&k); // k的范围为10**12 所以用long long读入
    for (int i = 1; i <= n; i++) { // 选择 index 从 1 开始
        for (int j = 1; j <= n; j++) {
            scanf("%d", &matrixA.coefficient[i][j]); // 初始矩阵上的系数均为int范围内, 因此%d读入不会爆int
        }
    }
    quickPower(k);
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            printf("%d ", result.coefficient[i][j]); // 结果经过mod所以%d不会爆int
        }
        fputs("\n", stdout);
    }
    return 0;
}