快速幂

定义

快速幂，二进制取幂（Binary Exponentiation，也称平方法），是一个在 $\mathrm{\Theta }(\log n)$ 的时间内计算 $a^n$ 的小技巧，而暴力的计算需要 $\mathrm{\Theta }(n)$ 的时间。

这个技巧也常常用在非计算的场景，因为它可以应用在任何具有结合律的运算中。其中显然的是它可以应用于模意义下取幂、矩阵幂等运算，我们接下来会讨论。

解释

计算 $a$ 的 $n$ 次方表示将 $n$ 个 $a$ 乘在一起：$a^n=\underset{n\text{ 个 a}}{\underbrace{a\times a\cdots \times a}}$。然而当 $a,n$ 太大的时侯，这种方法就不太适用了。不过我们知道：$a^{b+c}=a^b\cdot a^c,a^{2b}=a^b\cdot a^b=(a^b{)}^2$。二进制取幂的想法是，我们将取幂的任务按照指数的 二进制表示 来分割成更小的任务。

过程

迭代版本

首先我们将 $n$ 表示为 2 进制，举一个例子：

$$3^{13}=3^{(1101{)}_2}=3^8\cdot 3^4\cdot 3^1$$

因为 $n$ 有 $\lfloor {\log }_{2}n\rfloor +1$ 个二进制位，因此当我们知道了 $a^1,a^2,a^4,a^8,\dots ,a^{2^{\lfloor {\log }_{2}n\rfloor }}$ 后，我们只用计算 $\mathrm{\Theta }(\log n)$ 次乘法就可以计算出 $a^n$。

于是我们只需要知道一个快速的方法来计算上述 3 的 $2^k$ 次幂的序列。这个问题很简单，因为序列中（除第一个）任意一个元素就是其前一个元素的平方。举一个例子：

$$\begin{array}{rl}{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{array}$$

因此为了计算 $3^{13}$，我们只需要将对应二进制位为 1 的整系数幂乘起来就行了：

$$3^{13}=6561\cdot 81\cdot 3=1594323$$

将上述过程说得形式化一些，如果把 $n$ 写作二进制为 $(n_t{n}_{t-1}\cdots n_1{n}_0{)}_2$，那么有：

$$n=n_t{2}^t+n_{t-1}{2}^{t-1}+n_{t-2}{2}^{t-2}+\cdots +n_1{2}^1+n_0{2}^0$$

其中 $n_i\in \{0,1\}$。那么就有

$$\begin{array}{rl}{a}^n& =(a^{n_t{2}^t+\cdots +n_0{2}^0})\\ \\ & =a^{n_0{2}^0}\times a^{n_1{2}^1}\times \cdots \times a^{n_t{2}^t}\end{array}$$

根据上式我们发现，原问题被我们转化成了形式相同的子问题的乘积，并且我们可以在常数时间内从 $2^i$ 项推出 $2^{i+1}$ 项。

这个算法的复杂度是 $\mathrm{\Theta }(\log n)$ 的，我们计算了 $\mathrm{\Theta }(\log n)$ 个 $2^k$ 次幂的数，然后花费 $\mathrm{\Theta }(\log n)$ 的时间选择二进制为 1 对应的幂来相乘。

递归版本

上述迭代版本中，由于 $2^{i+1}$ 项依赖于 $2^i$，使得其转换为递归版本比较困难（一方面需要返回一个额外的 $a^{2^i}$，对函数来说无法实现一个只返回计算结果的接口；另一方面则是必须从低位往高位计算，即从高位往低位调用，这也造成了递归实现的困扰），下面则提供递归版本的思路。

给定形式 $n_{t\dots i}=(n_t{n}_{t-1}\cdots n_i{)}_2$，即 $n_{t\dots i}$ 表示将 $n$ 的前 $t-i+1$ 位二进制位当作一个二进制数，则有如下变换：

$$\begin{array}{rl}n& =n_{t\dots 0}\\ & =2\times n_{t\dots 1}+n_0\\ & =2\times (2\times n_{t\dots 2}+n_1)+n_0\\ & \cdots \end{array}$$

那么有：

$$\begin{array}{rl}{a}^n& =a^{n_{t\dots 0}}\\ & =a^{2\times n_{t\dots 1}+n_0}={\left(a^{n_{t\dots 1}}\right)}^2{a}^{n_0}\\ & ={\left(a^{2\times n_{t\dots 2}+n_1}\right)}^2{a}^{n_0}={\left({\left(a^{n_{t\dots 2}}\right)}^2{a}^{n_1}\right)}^2{a}^{n_0}\\ & \cdots \end{array}$$

如上所述，在递归时，对于不同的递归深度是相同的处理：$a^{n_{t\dots i}}=(a^{n_{t\dots (i+1)}}{)}^2{a}^{n_i}$，即将当前递归的二进制数拆成两部分：最低位在递归出来时乘上去，其余部分则变成新的二进制数递归进入更深一层作相同的处理。

可以观察到，每递归深入一层则二进制位减少一位，所以该算法的时间复杂度也为 $\mathrm{\Theta }(\log n)$。

应用

模意义下取幂

这是一个非常常见的应用，例如它可以用于计算模意义下的乘法逆元。既然我们知道取模的运算不会干涉乘法运算，因此我们只需要在计算的过程中取模即可。

首先我们可以直接按照上述递归方法实现：

1
2
3
4
5
6
7
8

long long binpow(long long a, long long b, long long p) {
  if (b == 0) return 1;
  long long res = binpow(a, b / 2, p);
  if (b % 2)
    return res * res % p * a % p;
  else
    return res * res % p;
}

1
2
3
4
5
6
7
8

def binpow(a, b, p):
    if b == 0:
        return 1
    res = binpow(a, b // 2, p)
    if (b % 2) == 1:
        return res * res * a % p
    else:
        return res * res % p

第二种实现方法是非递归式的。它在循环的过程中将二进制位为 1 时对应的幂累乘到答案中。尽管两者的理论复杂度是相同的，但第二种在实践过程中的速度是比第一种更快的，因为递归会花费一定的开销。

1
2
3
4
5
6
7
8
9

long long binpow(long long a, long long b, long long p) {
  long long res = 1;
  while (b > 0) {
    if (b & 1) res = res * a % p;
    a = a * a % p;
    b >>= 1;
  }
  return res;
}

1
2
3
4
5
6
7
8

def binpow(a, b, p):
    res = 1
    while b > 0:
        if b & 1:
            res = res * a % p
        a = a * a % p
        b >>= 1
    return res

计算斐波那契数

根据斐波那契数列的递推式 $F_n=F_{n-1}+F_{n-2}$，我们可以构建一个 $2\times 2$ 的矩阵来表示从 $F_i,F_{i+1}$ 到 $F_{i+1},F_{i+2}$ 的变换。于是在计算这个矩阵的 $n$ 次幂的时侯，我们使用快速幂的思想，可以在 $\mathrm{\Theta }(\log n)$ 的时间内计算出结果。对于更多的细节参见 斐波那契数列，矩阵快速幂的实现参见 矩阵加速递推 中的实现。

多次置换

简单地把这个置换取 $k$ 次幂，然后把它应用到序列上即可。时间复杂度为 $O(n\log k)$。对于更多的细节参见 置换的复合。

加速几何中对点集的操作

让我们来观察一下这三种操作对坐标的影响：

1. Shift 操作：将每一维的坐标分别加上一个常量；
2. Scale 操作：把每一维坐标分别乘上一个常量；
3. Rotate 操作：参考 Rotation matrix from axis and angle，可以用旋转矩阵乘以点对应的向量得到。

可以看到，每一个操作可以用一个 $4\times 4$ 的矩阵来表示：

$$\left[\begin{array}{cccc}{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{array}\right]$$

使用这个矩阵就可以将一个坐标（向量）进行变换，得到新的坐标（向量）：

$$\left[\begin{array}{cccc}{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{array}\right]\cdot \left[\begin{array}{c}x\\ y\\ z\\ 1\end{array}\right]=\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ z^{\prime }\\ 1\end{array}\right]$$

在矩阵中，我们在三维坐标之外添加了一维，并将其置为 $1$。这样的话，我们就可以使用矩阵乘法来描述 Shift 操作了。

对于每个操作，可以按如下方法构造矩阵。

1. （Shift 操作）让 $x$ 坐标方向的位移为 $\mathrm{\Delta }x$，$y$ 坐标的位移为 $\mathrm{\Delta }y$，$z$ 坐标的位移为 $\mathrm{\Delta }z$：

   $$\left[\begin{array}{cccc}1& 0& 0& \mathrm{\Delta }x\\ 0& 1& 0& \mathrm{\Delta }y\\ 0& 0& 1& \mathrm{\Delta }z\\ 0& 0& 0& 1\end{array}\right]$$

2. （Scale 操作）把 $x$ 坐标拉伸 $\mathrm{\Delta }x$ 倍，$y$ 坐标拉伸 $\mathrm{\Delta }y$ 倍，$z$ 坐标拉伸 $\mathrm{\Delta }z$ 倍：

   $$\left[\begin{array}{cccc}\mathrm{\Delta }x& 0& 0& 0\\ 0& \mathrm{\Delta }y& 0& 0\\ 0& 0& \mathrm{\Delta }z& 0\\ 0& 0& 0& 1\end{array}\right]$$

3. （Rotate 操作）绕单位向量 $\mathbf{u}=(u_x,u_y,u_z)$ 逆时针旋转 $\theta$ 弧度

   $$\left[\begin{array}{cccc}{u}_x^2(1-\cos \theta )+\cos \theta & u_x{u}_y(1-\cos \theta )-u_z\sin \theta & u_x{u}_z(1-\cos \theta )+u_y\sin \theta & 0\\ u_x{u}_y(1-\cos \theta )+u_z\sin \theta & u_y^2(1-\cos \theta )+\cos \theta & u_y{u}_z(1-\cos \theta )-u_x\sin \theta & 0\\ u_x{u}_z(1-\cos \theta )-u_y\sin \theta & u_y{u}_z(1-\cos \theta )+u_x\sin \theta & u_z^2(1-\cos \theta )+\cos \theta & 0\\ 0& 0& 0& 1\end{array}\right]$$

现在，每一种操作都被表示为了一个矩阵，变换序列可以用矩阵的乘积来表示，而一个 Repeat 操作相当于取一个矩阵的 $k$ 次幂。这样可以用 $O(m\log k)$ 的时间计算出整个变换序列最终形成的矩阵。最后将它应用到 $n$ 个点上，总复杂度 $O(n+m\log k)$。

定长路径计数

我们把该图的邻接矩阵 $M$ 取 $k$ 次幂，那么 $M_{i,j}$ 就表示从 $i$ 到 $j$ 长度为 $k$ 的路径的数目。该算法的复杂度是 $O(n^3\log k)$。有关该算法的细节请参见 矩阵 页面。

模意义下的整数乘法

与二进制取幂的思想一样，这次我们将其中的一个乘数表示为若干个 2 的整数次幂的和的形式。因为在对一个数做乘 2 并取模的运算的时侯，我们可以转化为加减操作防止整型溢出。这样可以在 $O(\log m)$ 的时间复杂度下解决问题。递归方法如下：

$$a\cdot b=\{\begin{array}{ll}0& \text{if }a=0\\ 2\cdot \frac{a}{2}\cdot b& \text{if }a>0\text{ and }a\text{ even}\\ 2\cdot \frac{a-1}{2}\cdot b+b& \text{if }a>0\text{ and }a\text{ odd}\end{array}$$

但在实际使用中，此方法由于引入了更大的计算复杂度导致时间效率不优。实际编程中通常利用 快速乘 来进行模数范围在 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
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

#include <cmath>
#include <cstring>
#include <iostream>
using namespace std;

const int M = 500;

int a[505], b[505], t[505];

// 大整数乘法
void mult(int x[], int y[]) {
  memset(t, 0, sizeof(t));
  for (int i = 1; i <= x[0]; i++) {
    for (int j = 1; j <= y[0]; j++) {
      if (i + j - 1 > M) 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 binpow(int p) {
  if (p == 1) {
    memcpy(b, a, sizeof(b));
    return;
  }
  binpow(p / 2);  // (2^(p/2))^2=2^p
  mult(b, b);     // 对 b 平方
  if (p % 2 == 1) mult(b, a);
}

int main() {
  cin.tie(nullptr)->sync_with_stdio(false);
  int p;
  cin >> p;
  a[0] = 1;  // 记录 a 数组的位数
  a[1] = 2;  // 对 2 进行平方
  b[0] = 1;  // 记录 b 数组的位数
  b[1] = 1;  // 答案数组
  binpow(p);
  cout << (int)(log10(2) * p) + 1 << '\n';
  b[1] -= 1;  // 最后一位减 1
  for (int i = M; i >= 1; i--) {
    cout << b[i];
    if ((i - 1) % 50 == 0) {
      cout << '\n';
    }
  }
}

同一底数与同一模数的预处理快速幂

在同一底数与同一模数的条件下，可以利用 分块思想，用一定的时间（一般是 $O(\sqrt{n})$）预处理后用 $O(1)$ 的时间回答一次幂询问。

过程

1. 选定一个数 $s$，预处理出 $a^0$ 到 $a^s$ 与 $a^{0\cdot s}$ 到 $a^{\lfloor \frac{p}{s}\rfloor \cdot s}$ 的值并存在两个数组里；
2. 对于每一次询问 $a^{b}modp$，将 $b$ 拆分成 $\lfloor {\displaystyle \frac{b}{s}}\rfloor \cdot s+bmods$，则 $a^b=a^{\lfloor \frac{b}{s}\rfloor \cdot s}\times a^{bmods}$，可以 $O(1)$ 求出答案。

关于这个数 $s$ 的选择，我们一般选择 $\sqrt{p}$ 或者一个大小适当的 $2$ 的次幂。选择 $\sqrt{p}$ 可以使预处理较优，选择 $2$ 的次幂可以使用位运算简化计算。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

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] % mod;
}

习题

• UVa 1230 - MODEX
• UVa 374 - Big Mod
• UVa 11029 - Leading and Trailing
• Codeforces - Parking Lot
• SPOJ - The last digit
• SPOJ - Locker
• SPOJ - Just add it

本页面部分内容译自博文 Бинарное возведение в степень 与其英文翻译版 Binary Exponentiation。其中俄文版版权协议为 Public Domain + Leave a Link；英文版版权协议为 CC-BY-SA 4.0。

本页面最近更新：2025/8/23 02:50:36，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：Ir1d, sshwy, CBW2007, Enter-tainer, Tiphereth-A, Xeonacid, ksyx, ouuan, Henry-ZHR, HeRaNO, iamtwz, luoguyuntianming, Marcythm, Peanut-Tang, Aquistcev, billchenchina, CCXXXI, chinggg, eyedeng, FFjet, Great-designer, H-J-Granger, hsfzLZH1, Hszzzx, Jude Gao, kenlig, kfy666, Konano, Menci, NachtgeistW, qwqAutomaton, shawlleyw, shenshuaijie, StudyingFather, TrisolarisHD, TRSWNCA, uqzjc, Zhoier
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用