快速幂

YZY Lv3

在算法的世界中,效率是衡量优劣的核心标尺。一个看似简单的操作,如求幂 ,其背后可能隐藏着巨大的优化空间。本文将系统性地剖析快速幂(Fast Exponentiation)算法,从其数学原理出发,深入探讨其在不同场景下的实现、变种及其与其他算法思想的深刻联系。我们的目标不仅是理解快速幂是什么,更是要掌握其思想精髓,并能灵活地将其应用到更广泛的问题域中。

1. 问题的源起:朴素的求幂

我们面临的第一个基本问题是:给定三个整数 ,计算 的值。

最直观的方法是模拟幂的定义,即 相乘。我们可以通过一个简单的循环来实现:

1
2
3
4
5
6
7
8
9
10
11
12
13
#include<bits/stdc++.h>
using namespace std;
using ll = long long;

// 题意概括:计算 (a^b) % m
// 朴素算法
ll naive_pow(ll a, ll b, ll m) {
ll res = 1;
for (int i = 0; i < b; ++i) {
res = res * a % m;
}
return res;
}
  • 时间复杂度: 。循环执行 次,每次进行一次乘法和取模操作。
  • 空间复杂度:

这种方法的简洁性是其唯一的优点。在 的值非常大时,例如 ,线性时间的复杂度是完全无法接受的。我们需要一种更高效的策略,能够将复杂度从线性级别显著降低。这便是快速幂算法登场的契机。

2. 核心思想:二进制分解

优化的关键在于打破“逐一相乘”的线性思维。我们能否一次性计算出更大步长的结果?答案是肯定的,其理论基石在于对指数 的二进制表示。

任何一个正整数 都可以唯一地表示为其二进制形式:

其中 在二进制表示下第 位的值。

根据指数定律 ,我们可以将 进行如下转换:

由于 只取 ,这个式子可以进一步简化。当 时,该项为 ;当 时,该项为 ,对乘积无影响。因此,我们只需要将那些 的二进制位为 的项 乘起来即可。

这个公式揭示了快速幂的核心:

  1. 底数的分解: 我们不再需要 ,而是需要一系列具有特定规律的底数:
  2. 递推关系: 这一系列底数可以通过反复平移得到:, , , …,
  3. 组合: 最终结果是这些 项的子集乘积,具体选择哪些项由指数 的二进制位决定。

让我们以计算 为例来直观感受这个过程。

  1. 指数 的二进制表示为
  2. 这意味着
  3. 因此,
  4. 我们来计算所需的各项:
  5. 的二进制位为 1 对应的项相乘:
    • 结果 =

整个计算过程的迭代次数取决于 的二进制位数,即 。这使得时间复杂度从 骤降至 ,实现了质的飞跃。

3. 算法实现与细节

基于二进制分解的思想,我们可以设计出两种等价的实现方式:递归和迭代。

3.1 递归实现

递归的思路非常自然地体现了分治思想。对于计算 ,我们有:

  • 基准情况: 若 ,则
  • 递归步骤:
    • 是偶数,则 。我们只需递归计算 ,然后将结果平方。
    • 是奇数,则 。由于 是偶数,问题转化为 。我们递归计算 ,平方后,再乘以一个

为了避免重复计算,我们将 的结果存储在一个临时变量中。

代码示例:递归快速幂

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#include<bits/stdc++.h>
using namespace std;
using ll = long long;

// 题意概括:计算 (a^b) % m
ll qmi_rec(ll a, ll b, ll m) {
if (b == 0) return 1 % m; // b=0时返回1,注意m=1的特殊情况
ll t = qmi_rec(a, b / 2, m); // 递归计算 a^(b/2)
t = t * t % m; // 平方
if (b % 2 == 1) t = t * a % m; // 如果b是奇数,多乘一个a
return t;
}

int main() {
ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
ll a, b, m;
cin >> a >> b >> m;
cout << qmi_rec(a, b, m) << endl;
return 0;
}
  • 时间复杂度: 。每次递归调用,指数 的规模大约减半。递归深度为
  • 空间复杂度: 。递归调用栈的深度。

3.2 迭代实现

尽管递归实现很优雅,但在竞赛中,迭代实现通常更受欢迎。它避免了函数调用的开销和潜在的栈溢出风险,并且空间复杂度为

迭代的实现直接模拟了我们在二进制分解示例中的手动过程。我们维护两个变量:一个结果 res (初始化为1),一个变化的底数 a。我们从低到高遍历 的二进制位:

  • 如果当前 的最低位是 1,说明最终结果的乘积中包含当前的底数项,于是 res = res * a
  • 无论当前位是 0 还是 1,底数都需要为下一轮做准备,即 a = a * a,从 变为
  • 然后我们将 右移一位(b >>= 1),相当于考察下一位。
  • 循环直到 变为 0。

代码示例:迭代快速幂

这是竞赛中最常用和最核心的快速幂模板。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#include<bits/stdc++.h>
using namespace std;
using ll = long long;

// 题意概括:计算 (a^b) % m
ll qmi(ll a, ll b, ll m) {
ll res = 1 % m; // 初始化结果,处理m=1的特殊情况
a %= m; // a可能大于m,先取模
while (b > 0) {
if (b & 1) res = res * a % m; // 若b的当前末位为1,则累乘
a = a * a % m; // 底数平方
b >>= 1; // b右移一位,考察下一位
}
return res;
}

int main() {
ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
ll a, b, m;
cin >> a >> b >> m;
cout << qmi(a, b, m) << endl;
return 0;
}
  • 时间复杂度: 。循环次数等于 的二进制位数。
  • 空间复杂度: 。只使用了常数个额外变量。

3.3 实现中的关键细节

  1. 取模运算 (Modulo Operation): 在计算 时,中间结果可能会非常大,超出标准整型甚至 long long 的表示范围。模运算的性质 是我们的护身符。为了防止溢出,每一次乘法运算后都必须立即进行取模。观察上面的代码,res = res * a % ma = a * a % m 都遵循了这一原则。
  2. 数据类型 (Data Types): 即使 都在 int 范围内,中间乘积 a * ares * a 仍可能溢出 int。例如,若 ,则 a*a 约为 ,超出了 int(约 )的范围。因此,在进行乘法运算时,相关变量(至少是 resa)必须声明为 long long 类型,以确保中间结果能够被安全存储。using ll = long long; 是一个非常好的实践。
  3. 初始值: res 初始化为 而非简单的 。这是为了处理一个极端但有效的测试用例:。任何数对 取模都应得到 ,但如果 res 初始化为 ,在循环未执行时(如 )会错误地返回 1 % m 能正确处理这种情况,当 时,结果为

3.4 O(log b) 乘法:快速乘

在进行模运算时,我们面临一个潜在的陷阱:即使 ,乘积 仍可能超出 long long 的表示范围(约 )。例如,当 都是 级别的数时,a * b 会溢出。为了解决这个问题,我们需要一种不直接计算 a * b 就能得到 的方法。

这个问题的解决方案,名为快速乘(或因其速度相对较慢而戏称为“龟速乘”),其思想与快速幂同源。我们将乘法 看作是 相加。

这与快速幂的二进制分解完全一致,只是将底层的“乘法”换成了“加法”,将单位元“1”换成了“0”。

代码示例:快速乘

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#include<bits/stdc++.h>
using namespace std;
using ll = long long;

// 题意概括: 计算 (a * b) % m,防止 a * b 溢出
ll mul(ll a, ll b, ll m) {
ll res = 0;
a %= m;
while (b > 0) {
if (b & 1) res = (res + a) % m; // 乘法变加法
a = (a + a) % m; // 底数平方变底数加倍
b >>= 1;
}
return res;
}
  • 时间复杂度:
  • 空间复杂度:

在 C++ 中,对于 64 位整数,我们还可以利用 __int128 类型来直接计算乘积,这在性能上通常优于 的快速乘。但快速乘的思想是独立于特定语言特性的,并且在处理更高精度或不同代数结构时,其原理依然适用。

1
2
3
4
// 使用 __int128 的版本,通常更快
ll mul_fast(ll a, ll b, ll m) {
return (ll)((__int128)a * b % m);
}

现在,我们可以将标准快速幂 qmi 函数中的 res * a % ma * a % m 替换为 mul(res, a, m)mul(a, a, m),从而写出一个能处理超大模数的快速幂版本。

4. 推广:矩阵快速幂

快速幂的威力远不止于数值计算。它的思想可以被抽象和推广。只要一种运算满足**结合律 (Associativity)**,我们就可以用快速幂来加速该运算的 次重复执行。

结合律定义为:,其中 代表某种运算。

矩阵乘法恰好满足结合律。这为我们解决一类特殊的线性递推关系问题打开了大门。

4.1 引入场景:斐波那契数列

斐波那契数列是线性递推的经典范例:

其中 。要求计算第 可能非常大),并对一个数 取模。

的递推算法在 很大时会超时。我们需要一个 的解法。这里的关键是找到状态转移的“常数”。

我们定义一个状态向量,它包含了计算下一状态所需的所有信息。对于斐波那契数列,状态向量可以是 。我们的目标是找到一个矩阵 ,使得:

根据递推公式 和恒等式 ,我们可以构造出这个转移矩阵 :

这个矩阵 就是斐波那契数列的“魔力常数”。

将此关系链式展开:

计算 的问题,现在被转化为了计算矩阵 次幂。

4.2 矩阵快速幂的实现

由于矩阵乘法满足结合律,我们可以直接套用快速幂的逻辑。我们需要:

  1. 定义矩阵结构: 一个 structclass 来表示一个 的矩阵。
  2. 实现矩阵乘法: 一个函数,输入两个矩阵,返回它们的乘积。注意,所有计算都要在模 的意义下进行。
  3. 实现矩阵快速幂: 一个函数,其结构与整数快速幂几乎完全相同,只是:
    • “乘法”操作被替换为矩阵乘法函数。
    • “单位元”(即 res 的初始值)从数字 1 变为单位矩阵(主对角线为1,其余为0的矩阵)。

代码示例:矩阵快速幂求斐波那契数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
#include<bits/stdc++.h>
using namespace std;
using ll = long long;

ll n, m;

// 矩阵结构体
struct Mat {
ll mat[2][2];
Mat() { // 构造函数,初始化为0
memset(mat, 0, sizeof mat);
}
};

// 矩阵乘法
Mat mul(Mat x, Mat y) {
Mat res;
for (int i = 0; i < 2; i++) {
for (int j = 0; j < 2; j++) {
for (int k = 0; k < 2; k++) {
res.mat[i][j] = (res.mat[i][j] + x.mat[i][k] * y.mat[k][j]) % m;
}
}
}
return res;
}

// 矩阵快速幂
Mat qmi(Mat a, ll b) {
Mat res;
// 初始化为单位矩阵
res.mat[0][0] = res.mat[1][1] = 1;

while (b > 0) {
if (b & 1) res = mul(res, a);
a = mul(a, a);
b >>= 1;
}
return res;
}

int main() {
ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
// 题意概括:计算斐波那契数列第n项模m的值 (n>=1)
// F_n = T^(n-1) * F_1
// F_{n-1} = ... * F_0
// 我们求F_n,所以需要 T^(n-1) 的第一行和 (F_1, F_0) 列向量的乘积
// F_1=1, F_0=0
cin >> n >> m;

if (n == 0) {
cout << 0 << endl;
return 0;
}
if (n == 1) {
cout << 1 % m << endl;
return 0;
}

Mat T;
T.mat[0][0] = 1; T.mat[0][1] = 1;
T.mat[1][0] = 1; T.mat[1][1] = 0;

Mat Tn = qmi(T, n - 1);

// F_n = Tn[0][0]*F_1 + Tn[0][1]*F_0 = Tn[0][0] * 1 + Tn[0][1] * 0
cout << Tn.mat[0][0] << endl;

return 0;
}
  • 时间复杂度: 。其中 是矩阵的维度。对于斐波那契数列, 是一个常数,所以复杂度为 来自于一次 矩阵乘法所需的三重循环。
  • 空间复杂度: 。存储矩阵所需的空间。

4.3 构造转移矩阵的一般方法

矩阵快速幂的难点通常不在于算法本身,而在于如何为给定的递推关系构造出正确的转移矩阵。
假设我们有一个 阶线性递推关系:

其中 是常数。

我们的状态向量需要包含 个连续项,才能递推到下一步。令状态向量为:

我们需要找到一个 的转移矩阵 使得
的形式是

我们逐行构造 :

  • 第一行: 对应 的计算。根据递推公式,。所以 的第一行就是
  • 第二行: 对应 的计算。这是一个恒等关系,。所以 的第二行是
  • 第三行: 对应 的计算。。所以 的第三行是
  • 第 k 行: 对应 的计算。。所以 的第 行是

最终得到的转移矩阵 如下:

有了这个通用的构造方法,我们就可以将任何常系数线性递推问题转化为矩阵快速幂来求解,只要递推的阶数 不太大。

5. 扩展技巧:处理带常数项的递推

我们已经掌握了处理形如 的齐次线性递推。但在实际问题中,递推式常常会包含一个常数项,例如 ,或者更复杂的

直接套用之前的矩阵构造方法会失效,因为常数项 游离于状态向量之外。解决办法是增广矩阵的思想:将常数项也纳入状态向量中。

考虑一个一般的 阶线性递推,带常数项

我们构造一个 维的状态向量,最后一维用来“携带”这个常数项:

相应的,我们需要一个 的转移矩阵 的形式是

我们来构造 使得

  • 第一行: 对应 的计算。。所以 的第一行是
  • 接下来 : 依然是恒等关系,如 。这部分与之前相同,只是矩阵维度增加了。第 行 () 是在第 列为1,其余为0。
  • 最后一行: 对应常数项 的“递推”。常数项是不变的,即 。所以 的最后一行是

这样构造出的转移矩阵 为:

有了这个增广后的转移矩阵,我们又回到了熟悉的世界:

其中 是初始状态向量。计算 即可用矩阵快速幂在 的时间内解决。

示例:求斐波那契数列前 n 项和

让我们看一个具体问题:计算 ,其中 是斐波那契数列。
我们有递推关系:

这里的 依赖于 ,而 本身也在递推。这提示我们需要一个包含 的状态向量。
令状态向量为
我们希望找到矩阵 使得

推导 的每一行:

  1. 。所以第一行是
  2. 。所以第二行是
  3. 。所以第三行是

转移矩阵
初始状态向量(当 时)为

我们只需用矩阵快速幂计算 ,然后与初始向量相乘即可得到

代码示例:斐波那契前缀和

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
#include<bits/stdc++.h>
using namespace std;
using ll = long long;

ll n, m;

struct Mat {
ll mat[3][3];
Mat() { memset(mat, 0, sizeof mat); }
};

Mat mul(Mat x, Mat y) {
Mat res;
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 3; j++) {
for (int k = 0; k < 3; k++) {
res.mat[i][j] = (res.mat[i][j] + x.mat[i][k] * y.mat[k][j]) % m;
}
}
}
return res;
}

Mat qmi(Mat a, ll b) {
Mat res;
res.mat[0][0] = res.mat[1][1] = res.mat[2][2] = 1; // 单位矩阵
while (b > 0) {
if (b & 1) res = mul(res, a);
a = mul(a, a);
b >>= 1;
}
return res;
}

int main() {
ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
// 题意概括: 求 S_n = sum_{i=1 to n} F_i mod m
cin >> n >> m;

if (n == 0) { cout << 0 << endl; return 0; }
if (n == 1) { cout << 1 % m << endl; return 0; }

Mat T;
T.mat[0][0] = 1; T.mat[0][1] = 1; T.mat[0][2] = 1;
T.mat[1][0] = 0; T.mat[1][1] = 1; T.mat[1][2] = 1;
T.mat[2][0] = 0; T.mat[2][1] = 1; T.mat[2][2] = 0;

Mat Tn = qmi(T, n - 1);

// V_n = T^(n-1) * V_1
// V_1 = (S_1, F_1, F_0)^T = (1, 1, 0)^T
// S_n = Tn[0][0]*S_1 + Tn[0][1]*F_1 + Tn[0][2]*F_0
// S_n = Tn[0][0]*1 + Tn[0][1]*1 + Tn[0][2]*0
ll ans = (Tn.mat[0][0] + Tn.mat[0][1]) % m;
cout << ans << endl;

return 0;
}
  • 时间复杂度: ,这里 是常数,所以是
  • 空间复杂度: ,即

这个技巧将矩阵快速幂的应用范围从齐次线性递推推广到了非齐次线性递推,极大地增强了它的实战能力。

6. 数论中的加速器:欧拉定理

快速幂本身是一种代数技巧,但当它与数论,特别是模算术结合时,其威力会得到进一步的放大。欧拉定理是这种协同作用的最佳体现。

6.1 费马小定理

费马小定理是欧拉定理在模数为素数时的特殊情况。它指出:
如果 是一个素数,且整数 不是 的倍数(即 ),则有:

这个定理的直接推论是,在模 的意义下,指数是以 为周期的。因此,对于任意正整数 ,我们有:

注意: 这个公式成立的前提是 为素数且 。当 时,,公式自然成立。当 的倍数时,,则 (对于 ),该公式可能不适用,需要特判。

6.2 欧拉定理

欧拉定理将费马小定理推广到了模数为任意正整数的情况。
首先,定义欧拉函数 :表示小于等于 的正整数中与 互质的数的个数。
例如,,因为与 10 互质的数有 1, 3, 7, 9。
是素数 ,则 都与 互质,所以

欧拉定理指出:
如果正整数 互质(即 ),则有:

这同样导出一个指数的简化规则:

这使得我们可以处理指数本身非常巨大的情况,比如 。我们可以先计算 ,得到一个较小的指数,然后再用快速幂求解。

6.3 扩展欧拉定理

欧拉定理要求 。当 不互质时,该定理不成立。然而,我们有更强的扩展欧拉定理:

这个分情况的公式看起来复杂,但它提供了一个普适的降幂法则。第三条是关键:当指数 足够大(不小于 )时,即使 不互质,指数部分的行为依然呈现出以 为周期的规律,只是需要加上一个 的偏移量来保证指数为正且大于等于

在竞赛中,为了编码方便,通常将三条合并为一个统一的形式来处理

  1. 计算指数
  2. 如果 ,则新指数为
  3. 如果 ,则新指数就是 本身。
  4. 最终计算

这个统一的过程完美覆盖了扩展欧拉定理的所有情况。

示例:

这是一个典型的应用场景,被称为“上帝的难题”的简化版。
题意概括: 给定 ,计算
可能非常大。

直接计算 是不可行的。根据扩展欧拉定理,我们只需要计算指数 取模后的值。
。我们需要计算
这又是一个快速幂问题!我们可以用 qmi(b, c, p) 来求得。
但是,根据扩展欧拉定理,我们需要判断 是否大于等于 。在 qmi(b, c, p) 的计算过程中,我们可以轻易地判断 的真实值是否超过了 。如果在计算中任何一步的中间结果已经大于等于 ,那么最终的 也一定大于等于

解题步骤:

  1. 计算 。这需要对 进行质因数分解。
  2. 设计一个特殊的快速幂函数 qmi_phi(b, c, p),用于计算 ,同时记录 的真实值是否大于等于
  3. 调用 qmi_phi 得到指数 exp 和一个标志位 flag
  4. 如果 flag 为真(表示 ),则最终的指数为 exp + p
  5. 否则,最终的指数就是 exp
  6. 用标准的快速幂 qmi(a, final_exp, m) 求解。

代码示例:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
#include<bits/stdc++.h>
using namespace std;
using ll = long long;

ll get_phi(ll n) {
ll res = n;
for (ll i = 2; i * i <= n; i++) {
if (n % i == 0) {
res = res / i * (i - 1);
while (n % i == 0) n /= i;
}
}
if (n > 1) res = res / n * (n - 1);
return res;
}

// flag用于标记b^c的真实值是否已经超过模数m
bool flag;
ll qmi_phi(ll a, ll b, ll m) {
ll res = 1;
flag = false;
a %= m;
while (b > 0) {
if (b & 1) {
if (res * a >= m) flag = true;
res = res * a % m;
}
if (b > 1 && a * a >= m) flag = true; // 提前判断,防止a*a溢出
a = a * a % m;
b >>= 1;
}
return res;
}

ll qmi(ll a, ll b, ll m) {
ll res = 1;
a %= m;
while (b > 0) {
if (b & 1) res = res * a % m;
a = a * a % m;
b >>= 1;
}
return res;
}

int main() {
ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
// 题意概括: 计算 a^(b^c) mod m
ll a, b, c, m;
cin >> a >> b >> c >> m;

ll p = get_phi(m);
ll exp = qmi_phi(b, c, p);
if (flag) { // 如果 b^c >= phi(m)
exp += p;
}

cout << qmi(a, exp, m) << endl;

return 0;
}
  • 时间复杂度: 。主要由计算 决定。
  • 空间复杂度:

7. 思想的延伸:二进制倍增

快速幂的核心思想——“二进制分解”和“倍增”——是一种极为深刻且普适的计算范式。它不仅限于代数运算,更可以应用于图论、数据结构等多个领域。其中,最著名的应用就是二进制倍增,尤其是在处理树上路径问题时,如求解最近公共祖先(LCA)。

让我们来剖析二者思想上的同构性:

快速幂 二进制倍增 (k-th Ancestor)
目标: 计算 次操作的累积结果。 目标: 从节点 向上跳 步。
操作: 乘法运算。 操作: 向上移动到父节点。
分解: 分解:
预处理/迭代: 计算 预处理: 计算 fa[u][i],即 的第 个祖先。
递推关系: 递推关系: fa[u][i] = fa[fa[u][i-1]][i-1]
组合: 若 的第 位为1,则结果乘以 组合: 若 的第 位为1,则从当前节点跳 步。

这个类比清晰地揭示了,快速幂求 的过程,本质上是从单位元 1 开始,根据 的二进制位,选择性地“乘以”预计算好的 次幂。
同样地,求 级祖先的过程,是从当前节点 开始,根据 的二进制位,选择性地“向上跳”预计算好的 步。

实现 k-th 祖先

  1. 预处理: 通过一次 DFS/BFS 得到每个节点的深度 dep[u] 和直接父节点 fa[u][0]

  2. 构建倍增表:

    1
    2
    3
    4
    5
    for (int j = 1; j < LOGN; ++j) {
    for (int i = 1; i <= n; ++i) {
    fa[i][j] = fa[fa[i][j-1]][j-1];
    }
    }

    这个双重循环的结构,完美地体现了 的思想。内层循环对所有节点 i 计算它们的 级祖先,这是通过它们的 级祖先的 级祖先得到的。

  3. 查询:

    1
    2
    3
    4
    5
    6
    // 寻找 u 的第 k 个祖先
    for (int i = LOGN - 1; i >= 0; --i) {
    if ((k >> i) & 1) { // 如果 k 的第 i 位是 1
    u = fa[u][i]; // 向上跳 2^i 步
    }
    }

    这个查询过程与迭代版快速幂的 while 循环异曲同工。

二进制倍增的思想还可以进一步扩展,用于在树上路径中维护具有可合并性的信息(如路径最大/最小值、路径和等)。这展示了快速幂所蕴含的“分治-合并”思想的强大生命力。

8. 变种与特殊应用

快速幂的思想不仅限于上述经典场景,通过巧妙的变形和嫁接,它能在更多领域发挥意想不到的作用。

8.1 空间换时间:预处理快速幂

场景: 存在大量查询,每次查询的底数 和模数 固定,只有指数 变化。
标准的快速幂对每次查询都进行 的计算。如果查询次数非常多(例如 次),总时间可能会过高。我们可以通过预处理,以空间换时间,将单次查询优化到接近

思想: 这是“分块”或“Meet-in-the-middle”思想的体现。我们将指数 进行分解。选择一个合适的块大小 (例如 ),将 表示为:

于是

  • 预处理:
    1. “小步” (Baby Steps): 计算并存储 。存入数组 lo
    2. “大步” (Giant Steps): 计算基底 。然后计算并存储 '_' allowed only in math mode(a_S)^0, (a_S)^1, \dots, (a_S)^{\lfloor \text{max_b}/S \rfloor} \pmod m。存入数组 hi
  • 查询: 对于给定的 ,计算 。结果就是 hi[q] * lo[r] % m

复杂度:

  • 预处理时间: '_' allowed only in math modeO(S + \text{max_b}/S)。当 '_' allowed only in math modeS \approx \sqrt{\text{max_b}} 时,预处理时间达到最优,约为 '_' allowed only in math modeO(\sqrt{\text{max_b}})
  • 查询时间:
  • 空间复杂度: '_' allowed only in math modeO(S + \text{max_b}/S),同样是 '_' allowed only in math modeO(\sqrt{\text{max_b}})

代码示例:预处理快速幂

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
#include<bits/stdc++.h>
using namespace std;
using ll = long long;

const int S = 1 << 16; // 块大小,65536
ll a, m;
vector<ll> lo, hi;

void pre() {
lo.resize(S);
hi.resize(S); // 假设 b 的最大值在 S*S 附近

lo[0] = 1;
for (int i = 1; i < S; ++i) {
lo[i] = lo[i-1] * a % m;
}

ll aS = 1;
// 使用快速幂计算 a^S
ll tmp = a;
int p = S;
while(p > 0) {
if(p & 1) aS = aS * tmp % m;
tmp = tmp * tmp % m;
p >>= 1;
}

hi[0] = 1;
for (int i = 1; i < S; ++i) {
hi[i] = hi[i-1] * aS % m;
}
}

ll query(ll b) {
ll q = b / S;
ll r = b % S;
return hi[q] * lo[r] % m;
}

int main() {
ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
// 题意概括: 多次查询 a^b % m,其中 a, m 固定
int T;
cin >> a >> m >> T;

pre(); // 预处理

while (T--) {
ll b;
cin >> b;
cout << query(b) << endl;
}
return 0;
}

8.2 超越内建类型:高精度快速幂

场景: 需要计算 的精确值,而结果的位数非常大,远远超过 long long 的表示范围。
此时,我们需要自己实现高精度算术,并将快速幂的框架应用于高精度数。

思想:

  1. 数据结构: 使用 std::vector<int>std::string 来存储大整数,每个元素代表一个数位或一个更大的“块”。
  2. 核心运算: 实现高精度乘法 big_mul(num1, num2)。这是整个算法的基石。
  3. 算法框架: 迭代快速幂的框架保持不变,只是其中的 1 变成了高精度数 “1”,乘法 * 变成了 big_mul

代码示例:高精度快速幂

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
#include<bits/stdc++.h>
using namespace std;
using ll = long long;

// 高精度数用string表示
string mul(string s1, string s2) {
int n1 = s1.size(), n2 = s2.size();
if (s1 == "0" || s2 == "0") return "0";
vector<int> res(n1 + n2, 0);
reverse(s1.begin(), s1.end());
reverse(s2.begin(), s2.end());

for (int i = 0; i < n1; i++) {
for (int j = 0; j < n2; j++) {
res[i + j] += (s1[i] - '0') * (s2[j] - '0');
}
}

for (int i = 0; i < n1 + n2 - 1; i++) {
if (res[i] >= 10) {
res[i+1] += res[i] / 10;
res[i] %= 10;
}
}

string ans;
int len = n1 + n2 - 1;
while (len > 0 && res[len] == 0) len--; // 去除前导零
for (int i = len; i >= 0; i--) {
ans += to_string(res[i]);
}
return ans;
}

string qmi(string a, int b) {
string res = "1";
while (b > 0) {
if (b & 1) res = mul(res, a);
a = mul(a, a);
b >>= 1;
}
return res;
}

int main() {
ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
// 题意概括: 计算大整数 a 的 b 次方
string a;
int b;
cin >> a >> b;
if (b == 0) {
cout << 1 << endl;
return 0;
}
cout << qmi(a, b) << endl;
return 0;
}
  • 时间复杂度: 设 len(N) 为大数 N 的长度。一次长度为 的大数乘法耗时 。在快速幂中,数的长度会快速增长,总复杂度分析较为复杂,但主要瓶颈在于 次高精度乘法。
  • 空间复杂度: 存储结果所需的空间,约为

8.3 几何变换的加速:矩阵与点集

场景: 对一个点集,重复施加同一组几何变换(如平移、旋转、缩放) 次,求所有点的最终位置。其中 可能非常大。

思想: 二维/三维的几何变换可以表示为矩阵乘法。这个转化的关键是引入**齐次坐标 (Homogeneous Coordinates)**。

  • 齐次坐标: 一个二维点 用三维列向量 表示。这个 “1” 是关键的增广维度。
  • 变换矩阵: 所有常见的仿射变换都可以表示为一个 3x3 矩阵。
    • 平移 :
    • 绕原点旋转 :
    • 绕原点缩放 :
  • 组合变换: 一系列变换 依次施加于点 上,等价于用组合矩阵 左乘点向量 。注意矩阵乘法的顺序是“从右到左”对应变换顺序。
  • 重复施加: 将该组合变换重复 次,等价于用矩阵 进行变换。

问题被完美地转化为了计算矩阵的 次幂,这正是矩阵快速幂的用武之地。

解题流程:

  1. 根据给定的几何变换序列,构造出每个变换的 3x3 矩阵(如果是三维变换则是 4x4)。
  2. 将这些矩阵按正确的顺序相乘,得到单次复合变换的总矩阵
  3. 使用矩阵快速幂计算最终的变换矩阵
  4. 对点集中的每个点,将其表示为齐次坐标向量 ,计算其最终位置

这种方法将一个对每个点都重复 次、总复杂度为 的过程,优化为了 的过程(其中 )。在 巨大时,实现了惊人的加速。

8.4 思想的抽象:幺半群上的快速幂

只要运算 满足结合律,就可以使用快速幂。一个集合 和一个其上的二元运算 如果满足结合律,则构成一个半群 (Semigroup)。如果这个半群还存在一个单位元 (Identity Element) (即对于任意 , ),则它构成一个**幺半群 (Monoid)**。

快速幂算法的本质,就是在幺半群上计算一个元素的 次幂。

  • res 的初始值是单位元
  • a = a * a 变成了 a = a \otimes a
  • res = res * a 变成了 res = res \otimes a

这提供了一个强大的视角,让我们能在更多抽象的场景中识别并应用快速幂。

应用实例:函数复合
考虑一个函数 。函数复合运算是结合律的:,所以 。幺半群是全体从 的函数集合,运算是函数复合,单位元是恒等函数
计算 次复合 ,即 ,就可以用快速幂来加速。我们不再是对数值进行幂运算,而是对“函数”这个实体。

这在解决一些涉及状态反复迭代的DP问题时特别有用。如果状态转移可以被描述为一个固定的线性变换(如矩阵),或者一个固定的函数,我们就可以用快速幂计算出“转移步”这个整体效果,从而实现 的状态跳转。

应用实例:DP优化
考虑一个DP问题,状态 只依赖于 ,且转移形式固定:,其中 是一个变换。
如果我们要从 计算 ,就需要应用 次变换

如果变换 的复合运算满足结合律,我们就可以用快速幂在 的时间内计算出 这个“终极变换”,然后一次性地应用到初始状态 上。矩阵快速幂优化线性递推就是这一思想最经典的体现。

8.5 逆问题:离散对数 (BSGS)

快速幂解决了 中的正向问题——已知 。而它的逆问题——已知 求最小的非负整数 ,则需要动用BSGS (Baby-Step Giant-Step) 算法。此算法与我们之前讨论的“预处理快速幂”思想一脉相承。

前提:
思想: 目标是寻找 使得 。我们再次使用分块思想,设 ,并将 表示为 ,其中

算法分为两步:

  1. Baby Steps: 计算所有 的值,其中 。将结果 {值: r} 存入一个哈希表(std::mapstd::unordered_map)。
  2. Giant Steps: 计算 。然后,依次计算 for from to 。对于每个计算出的值,在哈希表中查找是否存在。如果存在,假设查到的值为 ,对应的指数为 ,那么我们就找到了一个解

第一个找到的解即为最小非负解。

代码示例:BSGS

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#include<bits/stdc++.h>
using namespace std;
using ll = long long;

// 题意概括: 给定 a, y, m, 求最小非负整数 b 使得 a^b = y (mod m)
// 要求 gcd(a, m) == 1
ll bsgs(ll a, ll y, ll m) {
if (y % m == 1 % m && m != 1) return 0; // 特判 b=0

map<ll, int> hs;
int s = sqrt(m) + 1;

// Baby steps
ll t = y % m;
for (int r = 0; r < s; ++r) {
hs[t] = r;
t = t * a % m;
}

// a_S = a^s
ll as = 1;
for(int i = 0; i < s; ++i) as = as * a % m;
if (as == 0) return -1; // a^s=0, 只有当y=0时有解,但y=0不互质

// Giant steps
t = as;
for (int q = 1; q <= s; ++q) {
if (hs.count(t)) {
ll res = (ll)q * s - hs[t];
if(res >= 0) return res;
}
t = t * as % m;
}

return -1; // No solution
}
  • 时间复杂度: (使用 map) 或 (使用 unordered_map)。主要开销在于填充和查询哈希表。
  • 空间复杂度: ,用于存储哈希表。

BSGS算法是数论中一个优雅的工具,它将指数搜索问题的时间复杂度从朴素的 降至 ,是时空

9. 写在最后

从一个简单的求幂问题出发,我们踏上了一段精彩的算法优化之旅。

  • 我们从 的朴素循环,跃升至 快速幂,其核心在于对指数的二进制分解
  • 我们见证了这种思想的强大泛化能力,从数值运算推广到矩阵乘法,从而优雅地解决了各类线性递推问题,甚至包括带常数项的复杂情况。
  • 我们探索了它与数论的深度融合,利用欧拉定理及其扩展形式,攻克了指数本身是巨大幂次的难题,实现了指数级别的“降维打击”。
  • 我们最终揭示了其思想的本质——二进制倍增,并看到了它在图论(如 LCA)等看似无关领域中的精彩应用,理解了不同算法背后相通的智慧。

算法学习的精髓,正在于此——从一个具体问题出发,理解其核心思想,将其抽象、泛化,最终形成一套普适的思维工具,用以解决更广阔世界中的新问题。

附录:例题选讲

P1226 【模板】快速幂

题目描述

给你三个整数 ,求

输入格式

输入只有一行三个整数,分别代表

输出格式

输出一行一个字符串 a^b mod p=s,其中 分别为题目给定的值, 为运算结果。

输入输出样例 #1

输入 #1

1
2 10 9

输出 #1

1
2^10 mod 9=7

说明/提示

样例解释

数据规模与约定

对于 的数据,保证

CF630I Parking Lot

题目描述

停车场共有 个停车位。共有 种品牌的汽车,每种汽车的数量都远大于停车位的数量。

该公司首席执行官认为,如果停车场有 恰好 个连续汽车的品牌相同,则停车场会更漂亮。

给定n的值,问有多少的方案使停车场满足条件。

输入格式

一行,包含一个整数 )。

输出格式

输出一个整数,即总方案数。

输入输出样例 #1

输入 #1

1
3

输出 #1

1
24

P1045 [NOIP 2003 普及组] 麦森数

题目描述

形如 的素数称为麦森数,这时 一定也是个素数。但反过来不一定,即如果 是个素数, 不一定也是素数。到 1998 年底,人们已找到了 37 个麦森数。最大的一个是 ,它有 909526 位。麦森数有许多重要应用,它与完全数密切相关。

任务:输入 ,计算 的位数和最后 位数字(用十进制高精度数表示)

输入格式

文件中只包含一个整数

输出格式

第一行:十进制高精度数 的位数。

行:十进制高精度数 的最后 位数字。(每行输出 位,共输出 行,不足 位时高位补

不必验证 是否为素数。

输入输出样例 #1

输入 #1

1
1279

输出 #1

1
2
3
4
5
6
7
8
9
10
11
386
00000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000
00000000000000104079321946643990819252403273640855
38615262247266704805319112350403608059673360298012
23944173232418484242161395428100779138356624832346
49081399066056773207629241295093892203457731833496
61583550472959420547689811211693677147548478866962
50138443826029173234888531116082853841658502825560
46662248318909188018470682222031405210266984354887
32958028878050869736186900714720710555703168729087

说明/提示

【题目来源】

NOIP 2003 普及组第四题

P1226 【模板】快速幂

题目概括

计算 的值。

思路分析

本题是快速幂算法最直接、最纯粹的应用场景,是检验模板正确性的标准题目。问题要求计算 ,其中 的值较大,朴素的 做法显然无法通过。

我们直接套用在第 3.2 节中详细阐述的迭代快速幂算法即可。该算法的核心思想是利用指数 的二进制表示,将 次乘法操作优化为 次。

关键实现细节

  1. 数据类型:题目中 的范围小于 ,在 C++ 中可以用 long long 类型来存储以确保中间乘积 res * a 不会溢出。例如,当 时,resa 的最大值接近 ,其乘积约为 ,恰好在 long long 的表示范围内。
  2. 取模操作:在算法的每一步乘法后,都必须立即进行 % p 操作,以保证中间结果始终小于 ,防止溢出并确保最终结果的正确性。
  3. 输出格式:题目要求以特定字符串格式 "a^b mod p=s" 输出,需要在主函数中处理。

代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#include<bits/stdc++.h>
using namespace std;
using ll = long long;

// 题意概括: 计算 (a^b) % p
ll qmi(ll a, ll b, ll p) {
ll res = 1 % p;
a %= p;
while (b > 0) {
if (b & 1) res = res * a % p;
a = a * a % p;
b >>= 1;
}
return res;
}

int main() {
ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
ll a, b, p;
cin >> a >> b >> p;
ll s = qmi(a, b, p);
cout << a << "^" << b << " mod " << p << "=" << s << endl;
return 0;
}
  • 时间复杂度: 。循环的次数取决于指数 的二进制位数。
  • 空间复杂度:

CF630I Parking Lot

题目概括

在一个有 个车位的停车场,用 4 种品牌的车填满所有车位。要求恰好存在一个长度为 的连续段,其内的车辆品牌完全相同。求满足此条件的方案总数。

思路分析

这是一道将组合计数与快速幂相结合的题目。核心在于通过分类讨论,构造出计算方案数的数学公式,而公式的计算则依赖于快速幂。

我们分步来构建方案:

  1. 确定长度为 的同色块的位置和颜色

    • 该同色块的起始位置可以从 。共有 种选择。
    • 该同色块的颜色可以是 4 种品牌之一。有 4 种选择。
  2. 处理同色块的邻居
    为了保证“恰好”一个长度为 的同色块,该块的左右邻居(如果存在)的颜色不能与该块相同。

    我们对同色块的位置进行分类讨论:

    • Case 1: 块在最左端 (起始于位置 1)。
      块占据 [1, n]。其右邻居是位置 n+1,颜色不能和块相同,有 3 种选择。剩下的 个车位可以任意选择 4 种颜色。
      方案数:

    • Case 2: 块在最右端 (起始于位置 )。
      块占据 [n-1, 2n-2]。其左邻居是位置 n-2,有 3 种选择。剩下的 个车位可任意选择。
      方案数:

    • Case 3: 块在中间 (起始于 )。
      这样的起始位置共有 个。
      对于每一个这样的位置,块的左右都有邻居。左邻居有 3 种选择,右邻居也有 3 种选择,共 种。
      剩下的 个车位可任意选择。
      方案数:

  3. 整合方案
    总方案数是上述所有情况之和,再乘以块的颜色选择(4种)。

    化简该式:

    由于 的值不大(),最终结果会超出标准 int 但在 long long 范围内。计算 需要使用快速幂。

代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
#include<bits/stdc++.h>
using namespace std;
using ll = long long;

// 无模数的快速幂
ll qmi(ll a, ll b) {
ll res = 1;
while (b > 0) {
if (b & 1) res = res * a;
a = a * a;
b >>= 1;
}
return res;
}

int main() {
ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
ll n;
cin >> n;
// 公式为 (36n - 12) * 4^(n-4)
// 也可以写成 (9n - 3) * 4^(n-3)
ll ans = (9 * n - 3) * qmi(4, n - 3);
cout << ans << endl;
return 0;
}
  • 时间复杂度: 。主要耗时在于计算 的幂次。
  • 空间复杂度:

P1045 [NOIP 2003 普及组] 麦森数

题目概括

对于给定的整数 ,计算 的十进制位数,以及其末尾 500 位的值。

思路分析

本题是快速幂在两个不同维度上的综合应用:数值计算与高精度计算。

  1. 计算位数
    一个正整数 的十进制位数等于
    我们要计算 的位数。当 很大时, 非常接近,它们的位数是相同的。
    因此,位数 =
    利用对数性质 ,可得:
    位数 =
    我们可以使用 C++ 中 <cmath> 库的 log10 函数(它返回 double 类型)来计算 ,然后取整加一。

  2. 计算末尾 500 位
    求一个数的末尾 500 位,等价于求该数对 取模。
    我们需要计算 。这可以分解为两步:先计算 ,再减 1。
    模数 极其巨大,无法用任何内建整型存储。因此,这必须通过高精度算术来完成。

    我们将实现一个高精度版本的快速幂:

    • 表示: 用 stringvector<int> 存储大数。
    • 运算: 核心是实现高精度乘法 mul(big_num, big_num)
    • 框架: 快速幂的迭代框架保持不变,只是其中的乘法操作全部替换为我们实现的高精度乘法。

    计算 的过程如下:

    • 初始化结果 res 为高精度数 “1”。
    • 初始化底数 base 为高精度数 “2”。
    • 循环遍历 的二进制位。如果当前位是 1,则 res = mul(res, base);无论如何,base = mul(base, base)
    • 由于我们只关心末 500 位,可以在每次高精度乘法后,只保留结果的末 500 位,以防止数字过长导致 TLE。这相当于在整个高精度运算中隐式地对 取模。

    最后,对计算出的 的末 500 位执行高精度减 1,并处理好输出格式(补前导零,分行输出)。

代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
#include<bits/stdc++.h>
using namespace std;
using ll = long long;

const int L = 500; // 只保留末尾500位

// 高精度乘法 (string版)
string mul(string s1, string s2) {
int n1 = s1.size(), n2 = s2.size();
vector<int> res(n1 + n2, 0);
reverse(s1.begin(), s1.end());
reverse(s2.begin(), s2.end());

for (int i = 0; i < n1; i++) {
for (int j = 0; j < n2; j++) {
res[i + j] += (s1[i] - '0') * (s2[j] - '0');
}
}
for (int i = 0; i < n1 + n2 - 1; i++) {
if (res[i] >= 10) {
res[i + 1] += res[i] / 10;
res[i] %= 10;
}
}

string ans;
int len = n1 + n2 - 1;
while (len > 0 && res[len] == 0) len--;
for (int i = len; i >= 0; i--) ans += to_string(res[i]);

// 截断,只保留末尾L位
if (ans.size() > L) ans = ans.substr(ans.size() - L);
return ans;
}

// 高精度快速幂
string qmi(string a, int b) {
string res = "1";
while (b > 0) {
if (b & 1) res = mul(res, a);
a = mul(a, a);
b >>= 1;
}
return res;
}

// 高精度减1
void sub1(string& s) {
int n = s.size() - 1;
while(n >= 0) {
if (s[n] > '0') {
s[n]--;
break;
} else {
s[n] = '9';
n--;
}
}
}


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

// 1. 计算位数
int digits = floor(p * log10(2)) + 1;
cout << digits << endl;

// 2. 高精度快速幂计算 2^P 的末500位
string ans = qmi("2", p);

// 3. 高精度减1
sub1(ans);

// 4. 格式化输出
// 补足前导0
while (ans.size() < L) ans = "0" + ans;

for (int i = 0; i < L; ++i) {
cout << ans[i];
if ((i + 1) % 50 == 0) cout << endl;
}

return 0;
}
  • 时间复杂度: 次高精度乘法,每次乘法最多是 的复杂度。
  • 空间复杂度: 。存储高精度数。
  • Title: 快速幂
  • Author: YZY
  • Created at : 2025-07-05 01:51:37
  • Updated at : 2025-07-05 01:51:37
  • Link: https://blog.dtoj.team/2025/07/05/快速幂/
  • License: 三金家的Y同学