快速幂 与 矩阵快速幂

因为一些事情得重新拾起算法了

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

矩阵乘法:

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

AAmpm*p 的矩阵,BBpnp*n 的矩阵,称 mnm*n 的矩阵 CCAABB 的乘积,则 CC 的第 ii 行第 jj 列的元素为

cij=k=1paikbkj=ai1b1j+ai2b2j++aipbpjc_{ij} = \sum_{k=1}^p a_{ik}b_{kj} = a_{i1}b_{1j} + a_{i2}b_{2j} + \cdots + a_{ip}b_{pj}

mm 为行,pp 为列

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

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

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

A=[a11a12a1na21a22a2nam1am2amn]}nm.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++

#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列的元素为
// $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为列
// 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;
}