二进制求幂法（Binary Exponentiation）
二进制求幂法（也称为平方求幂）是一种高级技巧，它可以用
O
(
log
n
)
O(\log n)
O(logn) 的时间计算出
a
n
a^{n}
an （初级方法需要
O
(
n
)
O(n)
O(n) ）。
它在许多与计算无关的任务中也有重要应用，因为它可以与任何具有关联性属性的操作一起使用：
(
X
⋅
Y
)
⋅
Z
=
X
⋅
(
Y
⋅
Z
)
(X \cdot Y) \cdot Z = X \cdot (Y \cdot Z)
(X⋅Y)⋅Z=X⋅(Y⋅Z)
最明显的是，这适用于模块化乘法、矩阵的乘法和其他问题，我将在下面进行讨论。
算法
求
a
n
a^{n}
an 可以简单地表示为对
a
a
a 做
n
−
1
n-1
n−1 次乘法：
a
n
=
a
⋅
a
⋅
…
⋅
a
a^{n} = a \cdot a \cdot \ldots \cdot a
an=a⋅a⋅…⋅a 。但是这种方法对于较大的
a
a
a 和
n
n
n 没有实用性。
a
b
+
c
=
a
b
⋅
a
c
a^{b+c} = a^b \cdot a^c
ab+c=ab⋅ac 和
a
2
b
=
a
b
⋅
a
b
=
(
a
b
)
2
a^{2b} = a^b \cdot a^b = (a^b)^2
a2b=ab⋅ab=(ab)2
二进制求幂的想法是：使用指数的二进制表示来拆分工计算。
把
n
n
n 用二进制表示，例如：
3
13
=
3
110
1
2
=
3
8
⋅
3
4
⋅
3
1
3^{13} = 3^{1101_2} = 3^8 \cdot 3^4 \cdot 3^1
313=311012=38⋅34⋅31
数n由位二进制数组成，所以如果提前知道
a
1
,
a
2
,
a
4
,
a
8
,
…
,
a
⌊
log
n
⌋
a^1, a^2, a^4, a^8, \dots, a^{\lfloor \log n \rfloor}
a1,a2,a4,a8,…,a⌊logn⌋ 的结果，则只需要执行
O
(
log
n
)
O(\log n)
O(logn) 的乘法即可。
所以我们只需要知道一个快速的方法来计算这些。幸运的是这很容易，因为序列中的一个元素只是前一个元素的平方。
3
1
=
3
3
2
=
(
3
1
)
2
=
3
2
=
9
3
4
=
(
3
2
)
2
=
9
2
=
81
3
8
=
(
3
4
)
2
=
8
1
2
=
6561
\begin{aligned} 3^1 &= 3 \\ 3^2 &= \left(3^1\right)^2 = 3^2 = 9 \\ 3^4 &= \left(3^2\right)^2 = 9^2 = 81 \\ 3^8 &= \left(3^4\right)^2 = 81^2 = 6561 \end{aligned}
31323438=3=(31)2=32=9=(32)2=92=81=(34)2=812=6561
所以要得到最终的答案
3
13
3^{13}
313 ，只需要将其中的三个相乘（跳过
3
2
3^{2}
32 因为相应的位在
n
n
n 的二进制数中为
0
0
0 ）：
3
13
=
6561
⋅
81
⋅
3
=
1594323
3^{13} = 6561 \cdot 81 \cdot 3 = 1594323
313=6561⋅81⋅3=1594323
该算法的最终复杂度为
O
(
log
n
)
O(\log n)
O(logn)。
以下递归方法表达了相同的想法：
a
n
=
{
1
n
=
0
(
a
n
2
)
2
n
>
0
,
n
为偶数
(
a
n
−
1
2
)
2
⋅
a
n
>
0
,
n
为奇数
a^n = \begin{cases} 1 & n = 0 \\ \left(a^{\frac{n}{2}}\right)^2 & n > 0 \text{ , } n \text{ 为偶数}\\ \left(a^{\frac{n - 1}{2}}\right)^2 \cdot a & n > 0 \text{ , } n \text{ 为奇数}\\ \end{cases}
an=⎩⎪⎪⎨⎪⎪⎧1(a2n)2(a2n−1)2⋅an=0n>0 , n 为偶数n>0 , n 为奇数
实现
首先是递归方法，它是递归公式的直接转化：
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;
}
第二种方法无需递归即可完成相同的任务。它计算循环中的所有幂，并将这些幂与相应的设置位相乘。尽管两种方法的复杂度相同，但这种方法在实践中会更快，因为我们没有递归调用的开销。
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;
}
应用
1.大指数模数的高效计算（Effective computation of large exponents modulo a number）
问题： 计算
x
n
m
o
d
m
x^n \bmod m
xnmodm 。这是一个很常见的操作。例如，它用于计算模乘逆（modular multiplicative inverse）。
解决方案： 因为我们知道模块运算符不会干扰乘法（
a
⋅
b
≡
(
a
m
o
d
m
)
⋅
(
b
m
o
d
m
)
(
m
o
d
m
)
a \cdot b \equiv (a \bmod m) \cdot (b \bmod m) \pmod m
a⋅b≡(amodm)⋅(bmodm)(modm))，我们可以直接使用相同的代码，只需将每个乘法替换为模乘：
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;
}
注意： 如果
m
m
m 是一个素数，我们可以通过计算
x
n
m
o
d
(
m
−
1
)
x^{n \bmod (m-1)}
xnmod(m−1) 代替
x
n
x ^ n
xn 来加快这个算法的速度。这严格遵循费马小定理。
2.斐波那契数的高效计算（Effective computation of Fibonacci numbers）
问题： 计算第
n
n
n 个斐波那契数。
解决方案： 有关更多详细信息，请参阅斐波那契数那篇文章。这里将只对算法进行概述。要计算下一个斐波那契数，只需要前两个数，如
F
n
=
F
n
−
1
+
F
n
−
2
F_n = F_{n-1} + F_{n-2}
Fn=Fn−1+Fn−2 。我们可以建立一个
2
×
2
2 \times 2
2×2 的矩阵描述这种转换：从
F
i
F_i
Fi 和
F
i
+
1
F_{i+1}
Fi+1 到
F
i
+
1
F_{i+1}
Fi+1 和
F
i
+
2
F_{i+2}
Fi+2 。例如，将此转换应用于
F
0
F_0
F0 和
F
1
F_1
F1 则会把它改成
F
1
F_1
F1 和
F
2
F_2
F2 。因此，我们可以将这个变换矩阵扩大为
n
n
n 阶，则寻找
F
n
F_n
Fn 的时间复杂度就为
O
(
log
n
)
O(\log n)
O(logn) 。
3.应用排列
k
k
k 次（Applying a permutation
k
k
k times）
问题： 给定一个长度为
n
n
n 的序列。对其应用给定的排列
k
k
k 次。
解决方案： 只需将排列提高到
k
k
k 次幂使用二进制求幂，然后将其应用于序列。时间复杂度为
O
(
n
log
k
)
O(n \log k)
O(nlogk) 。
注意： 通过构建排列图并独立考虑每个循环，可以在线性时间内更有效地解决此任务。然后你可以计算
k
k
k 和循环大小的模，并找到作为该循环一部分的每个数字的最终位置。
4.在一组点上快速应用一组几何操作（Fast application of a set of geometric operations to a set of points）
问题： 给出
n
n
n 个点
p
i
p_i
pi ，对其中每个点应用
m
m
m 个变换。每个变换都可以是移位、缩放或旋转一个给定的角度。还有一个 "循环 "操作，它应用一个给定的变换列表
k
k
k 次（"循环 "操作可以嵌套）。你应该以快于
O
(
n
⋅
l
e
n
g
t
h
)
O(n \cdot length)
O(n⋅length) 的速度应用所有的变换，其中
l
e
n
g
t
h
length
length 是要应用的变换总数（结束 "循环 "操作后）。
解决方案： 让我们看看不同类型的变换是如何改变坐标的：
移位操作：为每个坐标添加不同的常数。缩放操作：将每个坐标乘以不同的常数。旋转操作：变换比较复杂（这里不再赘述），但是每个新坐标仍然可以表示为旧坐标的线性组合。
如您所见，每个变换都可以表示为坐标上的线性运算。因此，变换可以写成
4
×
4
4 \times 4
4×4 的矩阵：
(
a
11
a
12
a
13
a
14
a
21
a
22
a
23
a
24
a
31
a
32
a
33
a
34
a
41
a
42
a
43
a
44
)
\begin{pmatrix} a_{11} & a_ {12} & a_ {13} & a_ {14} \\ a_{21} & a_ {22} & a_ {23} & a_ {24} \\ a_{31} & a_ {32} & a_ {33} & a_ {34} \\ a_{41} & a_ {42} & a_ {43} & a_ {44} \end{pmatrix}
⎝⎜⎜⎛a11a21a31a41a12a22a32a42a13a23a33a43a14a24a34a44⎠⎟⎟⎞
即，当与一个具有旧坐标和单位的向量相乘时，会得到一个具有新坐标和单位的新向量：
(
x
y
z
1
)
⋅
(
a
11
a
12
a
13
a
14
a
21
a
22
a
23
a
24
a
31
a
32
a
33
a
34
a
41
a
42
a
43
a
44
)
=
(
x
′
y
′
z
′
1
)
\begin{pmatrix} x & y & z & 1 \end{pmatrix} \cdot \begin{pmatrix} a_{11} & a_ {12} & a_ {13} & a_ {14} \\ a_{21} & a_ {22} & a_ {23} & a_ {24} \\ a_{31} & a_ {32} & a_ {33} & a_ {34} \\ a_{41} & a_ {42} & a_ {43} & a_ {44} \end{pmatrix} = \begin{pmatrix} x' & y' & z' & 1 \end{pmatrix}
(xyz1)⋅⎝⎜⎜⎛a11a21a31a41a12a22a32a42a13a23a33a43a14a24a34a44⎠⎟⎟⎞=(x′y′z′1)
(你问，为什么要引入一个虚构的第四个坐标？没有这个，就不可能实现移位操作，因为它要求我们在坐标上添加一个常数。如果没有虚构的坐标，我们就只能对坐标进行线性组合，而不能增加一个常数。）
以下是一些如何以矩阵形式表示转换的示例：
移位操作：将
x
x
x 坐标移动
5
5
5 ，
y
y
y 坐标移动
7
7
7 ，
z
z
z 坐标移动
9
9
9 。
(
1
0
0
0
0
1
0
0
0
0
1
0
5
7
9
1
)
\begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 5 & 7 & 9 & 1 \end{pmatrix}
⎝⎜⎜⎛1005010700190001⎠⎟⎟⎞
缩放操作：
X
X
X 坐标的比例为
10
10
10 ，其他两个比例为
5
5
5 。
(
10
0
0
0
0
5
0
0
0
0
5
0
0
0
0
1
)
\begin{pmatrix} 10 & 0 & 0 & 0 \\ 0 & 5 & 0 & 0 \\ 0 & 0 & 5 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}
⎝⎜⎜⎛10000050000500001⎠⎟⎟⎞
旋转操作：按照右手规则（逆时针方向）围绕
X
X
X 轴旋转
θ
θ
θ 度。
(
1
0
0
0
0
cos
θ
−
sin
θ
0
0
sin
θ
cos
θ
0
0
0
0
1
)
\begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & \cos \theta & -\sin \theta & 0 \\ 0 & \sin \theta & \cos \theta & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}
⎝⎜⎜⎛10000cosθsinθ00−sinθcosθ00001⎠⎟⎟⎞
现在，一旦每个变换被描述为一个矩阵，变换的序列就可以被描述为这些矩阵的乘积，而一个重复了
k
k
k 次的 "循环 "可以被描述为矩阵提高到
k
k
k 的幂（这可以用
O
(
log
k
)
O(\log{k})
O(logk) 的二进制指数计算）。这样，代表所有变换的矩阵可以先在
O
(
m
log
k
)
O(m \log{k})
O(mlogk) 内计算出来，然后可以在
O
(
n
)
O(n)
O(n) 内应用于
n
n
n 个点中的每个点，总的复杂度为
O
(
n
+
m
log
k
)
O(n + m \log{k})
O(n+mlogk) 。
5.图中长度为
k
k
k 的路径数（Number of paths of length
k
k
k in a graph）
问题： 给定一个有
n
n
n 个顶点的无权图，找出从任何顶点
u
u
u 到任何其他顶点
v
v
v 的长度为
k
k
k 的路径数。
解决方案： 在另一篇文章中更详细地考虑了这个问题。该算法包括将图的邻接矩阵
M
M
M （一个矩阵，如果有一条从
i
i
i 到
j
j
j 的边，则
m
i
j
=
1
m_{ij}=1
mij=1 ，否则为
0
0
0 ）提高到
k
k
k 次方。现在，
m
i
j
m_{ij}
mij 将是长度为
k
k
k 的从
i
i
i 到
j
j
j 的路径数。该解决方案的时间复杂度为
O
(
n
3
log
k
)
O(n^3 \log k)
O(n3logk) 。
注意： 在同一篇文章中，考虑了这个问题的另一种变体：当边缘被加权时，需要找到恰好包含的最小权重路径边缘。如那篇文章所示，这个问题也可以通过邻接矩阵的取幂来解决。 矩阵将有从
i
i
i 到
j
j
j 的边的权重，如果没有这样的边，则为
∞
\infty
∞ 。代替通常的两个矩阵相乘的操作，应该使用一个修改过的：代替乘法，两个值相加，而不是求和，取最小值。结果为：
r
e
s
u
l
t
i
j
=
min
1
≤
k
≤
n
(
a
i
k
+
b
k
j
)
result_{ij} = \min\limits_{1\ \leq\ k\ \leq\ n}(a_{ik} + b_{kj})
resultij=1 ≤ k ≤ nmin(aik+bkj) 。
6.二进制求幂的变体：两个数字乘以
m
m
m 的模数（Variation of binary exponentiation: multiplying two numbers modulo
m
m
m ）
问题：
a
a
a 和
b
b
b 适合内置数据类型，但它们的乘积太大，不能放在64位的整数中，两个数字
a
a
a 和
b
b
b 乘以
m
m
m 。我们的想法是在不求较大值的情况下计算
a
⋅
b
(
m
o
d
m
)
a \cdot b \pmod m
a⋅b(modm) 。
解决方案： 我们简单地应用上面描述的二进制构造算法，只执行加法而不是乘法。换句话说，我们将两个数的乘法“扩展”为
O
(
log
m
)
O (\log m)
O(logm) 的加法和乘以二的运算（本质上是加法）。
a
⋅
b
=
{
0
a
=
0
2
⋅
a
2
⋅
b
a
>
0
,
a
为偶数
2
⋅
a
−
1
2
⋅
b
+
b
a
>
0
,
a
为奇数
a \cdot b = \begin{cases} 0 &a = 0 \\ 2 \cdot \frac{a}{2} \cdot b &a > 0 \text{ , }a \text{为偶数} \\ 2 \cdot \frac{a-1}{2} \cdot b + b &a > 0 \text{ , }a \text{为奇数} \end{cases}
a⋅b=⎩⎪⎨⎪⎧02⋅2a⋅b2⋅2a−1⋅b+ba=0a>0 , a为偶数a>0 , a为奇数
注意： 你可以通过使用浮点运算，以不同的方式解决这个任务。首先用浮点数计算表达式
a
⋅
b
m
\frac{a \cdot b}{m}
ma⋅b ，并将其转换为无符号整数
q
q
q ，用无符号整数运算法则从
a
⋅
b
a\cdot b
a⋅b 中减去
q
⋅
m
q\cdot m
q⋅m ，然后取
m
m
m 的模，就可以找到答案。这个解决方案看起来相当不可靠，但它非常快，而且非常容易实现。更多信息请看这里。