快速幂
快速幂,二进制取幂(Binary Exponentiation,也称平方法),是一个在 的时间内计算 的小技巧,而暴力的计算需要 的时间。而这个技巧也常常用在非计算的场景,因为它可以应用在任何具有结合律的运算中。其中显然的是它可以应用于模意义下取幂、矩阵幂等运算,我们接下来会讨论。
算法描述
计算 的 次方表示将 个 乘在一起:。然而当 太大的时侯,这种方法就不太适用了。不过我们知道:。二进制取幂的想法是,我们将取幂的任务按照指数的 二进制表示 来分割成更小的任务。
首先我们将 表示为 2 进制,举一个例子:
因为 有 个二进制位,因此当我们知道了 后,我们只用计算 次乘法就可以计算出 。
于是我们只需要知道一个快速的方法来计算上述 3 的 次幂的序列。这个问题很简单,因为序列中(除第一个)任意一个元素就是其前一个元素的平方。举一个例子:
因此为了计算 ,我们只需要将对应二进制位为 1 的整系数幂乘起来就行了:
将上述过程说得形式化一些,如果把 写作二进制为 ,那么有:
其中 。那么就有
根据上式我们发现,原问题被我们转化成了形式相同的子问题的乘积,并且我们可以在常数时间内从 项推出 项。
这个算法的复杂度是 的,我们计算了 个 次幂的数,然后花费 的时间选择二进制为 1 对应的幂来相乘。
代码实现
首先我们可以直接按照上述递归方法实现:
| // C++ Version
long long binpow(long long a, long long b) {
if (b == 0) return 1;
long long res = binpow(a, b / 2);
if (b % 2)
return res * res * a;
else
return res * res;
}
|
| # Python Version
def binpow(a, b):
if b == 0:
return 1
res = binpow(a, b // 2)
if (b % 2) == 1:
return res * res * a
else:
return res * res
|
第二种实现方法是非递归式的。它在循环的过程中将二进制位为 1 时对应的幂累乘到答案中。尽管两者的理论复杂度是相同的,但第二种在实践过程中的速度是比第一种更快的,因为递归会花费一定的开销。
| // C++ Version
long long binpow(long long a, long long b) {
long long res = 1;
while (b > 0) {
if (b & 1) res = res * a;
a = a * a;
b >>= 1;
}
return res;
}
|
| # Python Version
def binpow(a, b):
res = 1
while b > 0:
if (b & 1):
res = res * a
a = a * a
b >>= 1
return res
|
模板:Luogu P1226
应用
模意义下取幂
问题描述
计算 。
这是一个非常常见的应用,例如它可以用于计算模意义下的乘法逆元。
既然我们知道取模的运算不会干涉乘法运算,因此我们只需要在计算的过程中取模即可。
| // C++ Version
long long binpow(long long a, long long b, long long m) {
a %= m;
long long res = 1;
while (b > 0) {
if (b & 1) res = res * a % m;
a = a * a % m;
b >>= 1;
}
return res;
}
|
| # Python Version
def binpow(a, b, m):
a = a % m
res = 1
while b > 0:
if (b & 1):
res = res * a % m
a = a * a % m
b >>= 1
return res
|
注意:根据费马小定理,如果 是一个质数,我们可以计算 来加速算法过程。
计算斐波那契数
问题描述
计算斐波那契数列第 项 。
根据斐波那契数列的递推式 ,我们可以构建一个 的矩阵来表示从 到 的变换。于是在计算这个矩阵的 次幂的时侯,我们使用快速幂的思想,可以在 的时间内计算出结果。对于更多的细节参见 斐波那契数列。
多次置换
问题描述
给你一个长度为 的序列和一个置换,把这个序列置换 次。
简单地把这个置换取 次幂,然后把它应用到序列 上即可。时间复杂度是 的。
注意:给这个置换建图,然后在每一个环上分别做 次幂(事实上做一下 对环长取模的运算即可)可以取得更高效的算法,达到 的复杂度。
加速几何中对点集的操作
三维空间中, 个点 ,要求将 个操作都应用于这些点。包含 3 种操作:
- 沿某个向量移动点的位置(Shift)。
- 按比例缩放这个点的坐标(Scale)。
- 绕某个坐标轴旋转(Rotate)。
还有一个特殊的操作,就是将一个操作序列重复 次(Loop),这个序列中也可能有 Loop 操作(Loop 操作可以嵌套)。现在要求你在低于 的时间内将这些变换应用到这个 个点,其中 表示把所有的 Loop 操作展开后的操作序列的长度。
让我们来观察一下这三种操作对坐标的影响:
- Shift 操作:将每一维的坐标分别加上一个常量;
- Scale 操作:把每一维坐标分别乘上一个常量;
- Rotate 操作:这个有点复杂,我们不打算深入探究,不过我们仍然可以使用一个线性组合来表示新的坐标。
可以看到,每一个变换可以被表示为对坐标的线性运算,因此,一个变换可以用一个 的矩阵来表示:
使用这个矩阵就可以将一个坐标(向量)进行变换,得到新的坐标(向量):
你可能会问,为什么一个三维坐标会多一个 1 出来?原因在于,如果没有这个多出来的 1,我们没法使用矩阵的线性变换来描述 Shift 操作。
接下来举一些简单的例子来说明我们的思路:
-
Shift 操作:让 坐标方向的位移为 , 坐标的位移为 , 坐标的位移为 :
-
Scale 操作:把 坐标拉伸 10 倍, 坐标拉伸 5 倍:
-
Rotate 操作:绕 轴旋转 弧度,遵循右手定则(逆时针方向)
现在,每一种操作都被表示为了一个矩阵,变换序列可以用矩阵的乘积来表示,而一个 Loop 操作相当于取一个矩阵的 k 次幂。这样可以用 计算出整个变换序列最终形成的矩阵。最后将它应用到 个点上,总复杂度 。
定长路径计数
问题描述
给一个有向图(边权为 1),求任意两点 间从 到 ,长度为 的路径的条数。
我们把该图的邻接矩阵 M 取 k 次幂,那么 就表示从 到 长度为 的路径的数目。该算法的复杂度是 。有关该算法的细节请参见 矩阵 页面。
模意义下大整数乘法
计算 。
与二进制取幂的思想一样,这次我们将其中的一个乘数表示为若干个 2 的整数次幂的和的形式。因为在对一个数做乘 2 并取模的运算的时侯,我们可以转化为加减操作防止溢出。这样仍可以在 的时内解决问题。递归方法如下:
快速乘
但是 的“龟速乘”还是太慢了,这在很多对常数要求比较高的算法比如 Miller_Rabin 和 Pollard-Rho 中,就显得不够用了。所以我们要介绍一种可以处理模数在 long long
范围内、不需要使用黑科技 __int128
的、复杂度为 的“快速乘”。
我们发现:
我们巧妙运用 unsigned long long
的自然溢出:
于是在算出 后,两边的乘法和中间的减法部分都可以使用 unsigned long long
直接计算,现在我们只需要解决如何计算 。
我们考虑先使用 long double
算出 再乘上 。
既然使用了 long double
,就无疑会有精度误差。极端情况就是第一个有效数字(二进制下)在小数点后一位。在 x86-64
机器下,long double
将被解释成 位拓展小数(即符号为 位,指数为 位,尾数为 位),所以 long double
最多能精确表示的有效位数为 。所以 最差从第 位开始出错,误差范围为 。乘上 这个 位整数,误差范围为 ,再加上 误差范围为 ,取整后误差范围位 。于是乘上 后,误差范围变成 ,我们需要判断这两种情况。
因为 在 long long
范围内,所以如果计算结果 在 时,直接返回 ,否则返回 ,当然你也可以直接返回 。
代码实现如下:
| long long binmul(long long a, long long b, long long m) {
unsigned long long c =
(unsigned long long)a * b -
(unsigned long long)((long double)a / m * b + 0.5L) * m;
if (c < m) return c;
return c + m;
}
|
高精度快速幂
前置技能
请先学习 高精度
例题【NOIP2003 普及组改编·麦森数】(原题在此)
题目大意:从文件中输入 P(1000<P<3100000),计算 的最后 100 位数字(用十进制高精度数表示),不足 100 位时高位补 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 | #include <bits/stdc++.h>
using namespace std;
int a[505], b[505], t[505], i, j;
int mult(int x[], int y[]) // 高精度乘法
{
memset(t, 0, sizeof(t));
for (i = 1; i <= x[0]; i++) {
for (j = 1; j <= y[0]; j++) {
if (i + j - 1 > 100) continue;
t[i + j - 1] += x[i] * y[j];
t[i + j] += t[i + j - 1] / 10;
t[i + j - 1] %= 10;
t[0] = i + j;
}
}
memcpy(b, t, sizeof(b));
}
void ksm(int p) // 快速幂
{
if (p == 1) {
memcpy(b, a, sizeof(b));
return;
}
ksm(p / 2); //(2^(p/2))^2=2^p
mult(b, b); // 对b平方
if (p % 2 == 1) mult(b, a);
}
int main() {
int p;
scanf("%d", &p);
a[0] = 1; // 记录a数组的位数
a[1] = 2; // 对2进行平方
b[0] = 1; // 记录b数组的位数
b[1] = 1; // 答案数组
ksm(p);
for (i = 100; i >= 1; i--) {
if (i == 1) {
printf("%d\n", b[i] - 1); // 最后一位减1
} else
printf("%d", b[i]);
}
}
|
同一底数与同一模数的预处理快速幂
在同一底数与同一模数的条件下,可以利用分块思想,用一定的时间(一般是 )预处理后用 的时间回答一次幂询问。
算法的具体步骤是:
- 选定一个数 ,预处理出 到 与 到 的值并存在一个(或两个)数组里;
- 对于每一次询问 ,将 拆分成 ,则 ,可以 求出答案。
关于这个数 的选择,我们一般选择 或者一个大小适当的 的次幂(选择 可以使预处理较优,选择 的次幂可以使用位运算优化/简化计算)。
参考代码
| int pow1[65536], pow2[65536];
void preproc(int a, int mod) {
pow1[0] = pow2[0] = 1;
for (int i = 1; i < 65536; i++) pow1[i] = 1LL * pow1[i - 1] * a % mod;
int pow65536 = 1LL * pow1[65535] * a % mod;
for (int i = 1; i < 65536; i++) pow2[i] = 1LL * pow2[i - 1] * pow65536 % mod;
}
int query(int pows) { return 1LL * pow1[pows & 65535] * pow2[pows >> 16]; }
|
习题
参考资料与注释
build本页面最近更新:,更新历史
edit发现错误?想一起完善? 在 GitHub 上编辑此页!
people本页面贡献者:OI-wiki
copyright本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供,附加条款亦可能应用