快速幂 与 矩阵快速幂
因为一些事情得重新拾起算法了
- Reference(Aka thanks)(可能"Reference"在这里的用法不太对,欢迎提醒我来纠正)
- Wiki
- 基础
- [[Binary Exponentiation]]
- [[Linear Algebra|矩阵乘法]]
- 视频
- 题目
- 讨论(想法)
- 似乎还可以通过剪枝优化,下次()
- 关于快速幂的名字笔者并无好感,更喜欢称其为“二进制求幂”,非常直观
把基础拿出来讲 -- 矩阵乘法
矩阵乘法:
TL;DR: 前一个矩阵逐行
与 后一个矩阵逐列
相乘求和 -- 结果放入结果矩阵 字排列
设 为 的矩阵, 为 的矩阵,称 的矩阵 为 与 的乘积,则 的第 行第 列的元素为
为行, 为列
- Complexity
- 朴素求幂复杂度为 ,根据主定理(<-朋友告诉的,笔者还不会主定理),分析得到快速幂的时间复杂度为 ,
实际上还可以通过费马小定理加速if isPrime(m)
- 进而来到矩阵快速幂的复杂度: 设矩阵维度为 ,矩阵乘法的复杂度为 ,快速幂的复杂度为 ,故矩阵快速幂的时间复杂度为 ,而 为常数,而根据 时间复杂度 - 维基百科,自由的百科全书 中所述,在大O表示法中常数不计。矩阵快速幂复杂度即为
- 朴素求幂复杂度为 ,根据主定理(<-朋友告诉的,笔者还不会主定理),分析得到快速幂的时间复杂度为 ,
初学编程的伙伴可能对于循环中 不是很敏感,或者容易弄混,这里特意排出这个矩阵,供参考:
一个 的矩阵是一个由 行 列元素排列成的矩形阵列。即形如
最开始题目所对应的答案(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;
}