06.快速幂

autor: iamtwz, billchenchina, CBW2007, CCXXXI, chinggg, Enter-tainer, eyedeng, FFjet, gaojude, Great-designer, H-J-Granger, Henry-ZHR, hsfzLZH1, Ir1d, kenlig, Konano, ksyx, luoguyuntianming, Marcythm, Menci, NachtgeistW, ouuan, Peanut-Tang, qwqAutomaton, sshwy, StudyingFather, Tiphereth-A, TrisolarisHD, TRSWNCA, Xeonacid, Yuuko10032, Zhangjiacheng2006, Zhoier, Hszzzx, shenshuaijie, kfy666

引入

快速幂（fast exponentiation），也称 二进制取幂（binary exponentiation）或 平方取幂法（exponentiation by squaring），是一个在 $Θ (lo g n)$ 的时间内计算 $a^{n}$ 的小技巧，而暴力的计算需要 $Θ (n)$ 的时间。

这个技巧可以应用于任何 $a$ 的乘法满足结合律的场景中，例如模意义下取幂、矩阵幂等，详见后文应用一节。

过程

计算 $a$ 的 $n$ 次方表示将 $n$ 个 $a$ 乘在一起： $a^{n} = n 个 a a \times a \dots \times a$ 。然而当 $n$ 太大或单次乘法开销太大的时侯，这种方法就不太适用了。二进制取幂的想法是，将取幂的任务按照指数的 二进制表示 来分割成更小的任务。

例子

假设要计算 $3^{13}$ 。如果将它展开为连乘式，需要 $13 - 1 = 12$ 次乘法。但是，因为
$3^{13} = 3^{(1101)_{2}} = 3^{8} \times 3^{4} \times 3^{1},$
所以，只要能快速计算出 $3^{1}, 3^{2}, 3^{4}, 3^{8}$ ，就能通过 $2$ 次乘法计算出 $3^{13}$ 的值。于是，只需要知道一个快速的方法来计算上述 $3$ 的 $2^{k}$ 次幂的序列。这是容易的，因为因为序列中（除第一个）任意一个元素都是其前一个元素的平方。

根据这些分析，可以得到 $3^{13}$ 的计算过程如下：
$3^{1} 3^{2} 3^{4} 3^{8} 3^{13} = 3, = (3^{1})^{2} = 3^{2} = 9, = (3^{2})^{2} = 9^{2} = 81, = (3^{4})^{2} = 8 1^{2} = 6561, = 6561 \times 81 \times 3 = 1594323.$
过程中，只进行了 $5$ 次乘法运算。

这就是快速幂的基本想法。至于具体实现，有两种常见的版本。

迭代版本

设 $n$ 的二进制表示为 $(n_{t} n_{t - 1} \dots n_{1} n_{0})_{2}$ ，也就是说，有

n = n_{t} 2^{t} + n_{t - 1} 2^{t - 1} + \dots + n_{1} 2^{1} + n_{0} 2^{0},

其中， $n_{i} \in {0, 1}$ 。那么，就有

a^{n} = a^{n_{t} 2^{t} + n_{t - 1} 2^{t - 1} + \dots + n_{1} 2^{1} + n_{0} 2^{0}} = a^{n_{0} 2^{0}} \times a^{n_{1} 2^{1}} \times \dots \times a^{n_{t - 1} 2^{t - 1}} \times a^{n_{t} 2^{t}} .

注意，只有 $n_{i} = 1$ 的项才会真正出现在乘积的计算中。

根据这一表达式，可以首先在 $Θ (lo g n)$ 时间内计算出 $a$ 的 $Θ (lo g n)$ 个 $2^{k}$ 次幂的取值，然后花费 $Θ (lo g n)$ 的时间选择等于 $1$ 的二进制位对应的幂次乘到最终结果中。这就是快速幂的迭代版本实现。

伪代码如下：

Algorithm FastPow (a, n) : Input. Base a and exponent n . Output. Power a^{n} . Method. 123456789 result \leftarrow Id while n > 0 do if n mod 2 = 1 then result \leftarrow result \cdot a end if a \leftarrow a \cdot a n \leftarrow n /2 end while return result

利用这一方法计算快速幂，需要进行 $Θ (lo g n)$ 次乘法运算。

递归版本

这一过程同样可以通过递归形式实现。注意到，指数 $n$ 的二进制展开可以递归地写作

(n_{t} n_{t - 1} \dots n_{1} n_{0})_{2} = 2 \times (n_{t} n_{t - 1} \dots n_{1})_{2} + n_{0} .

因此，幂次 $a^{n}$ 可以递归地计算为

a^{n} = ⎩ ⎨ ⎧ 1, (a^{⌊ n /2 ⌋})^{2}, (a^{⌊ n /2 ⌋})^{2} \cdot a, n = 0, n > 0 and n is even, n > 0 and n is odd .

这就是快速幂的递归版本实现。

伪代码如下：

Algorithm FastPow (a, n) : Input. Base a and exponent n . Output. Power a^{n} . Method. 123456789 if n = 0 then return Id end if result \leftarrow FastPow (a, n /2) if n mod 2 = 0 then return result \cdot result \cdot a else return result \cdot result end if

利用这一方法计算快速幂，需要递归 $Θ (lo g n)$ 次，同样需要 $Θ (lo g n)$ 次乘法运算。尽管复杂度相同，由于递归本身有一定开销，所以实践中迭代版本的速度更快。

应用

模意义下取幂

洛谷 P1226【模板】快速幂

给定三个整数 $a, b, p$ ，求 $a^{b} mod p$ 。其中 $p \geq 2$ 。

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

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

参考实现

[list2tab]

C++

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;
}

Python

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 时对应的幂累乘到答案中。尽管两者的理论复杂度是相同的，但第二种在实践过程中的速度是比第一种更快的，因为递归会花费一定的开销。

参考实现

[list2tab]

C++

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;
}

Python

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

[!warning] 注意

模数通常情况下大于 $1$ 。在十分特殊的情况下，模数 $p$ 可能等于 $1$ ，此时需要特殊考虑 $b = 0$ 的情况。
当指数很大时，需利用扩展欧拉定理降幂后计算。

计算斐波那契数

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

多次置换

问题描述

给你一个长度为 $n$ 的序列和一个置换，把这个序列置换 $k$ 次。

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

注意

对这个置换建图，然后在每一个环上分别做 $k$ 次幂（事实上等价于 $k$ 对环长取模），可以在 $O (n)$ 的时间复杂度下解决此问题。

加速几何中对点集的操作

HDU 4087 A Letter to Programmers

给定三维空间中 $n$ 个点 $p_{i}$ ，要求将 $m$ 个操作都应用于这些点。包含 3 种操作：

沿某个向量移动点的位置（Shift）。

按比例缩放这个点的坐标（Scale）。

绕某条直线旋转（Rotate）。

还有一个特殊的操作，就是将某个操作序列重复 $k$ 次（Repeat），Repeat 操作可以嵌套。输出操作结束后每个点的坐标。

参考向量与矩阵中的内容，每一种操作都可以用一个变换矩阵表示，一系列连续的变换可以用矩阵的乘积来表示。一个 Repeat 操作就相当于取一个矩阵的 $k$ 次幂。这样可以用 $O (m lo g k)$ 的时间计算出整个变换序列最终形成的矩阵。最后将它应用到 $n$ 个点上，总复杂度 $O (n + m lo g k)$ 。

定长路径计数

问题描述

给一个有向图（边权为 1），求任意两点 $u, v$ 间从 $u$ 到 $v$ ，长度为 $k$ 的路径的条数。

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

模意义下的整数乘法

问题描述

给定非负整数 $a, b$ 和正整数 $m$ ，计算 $a \times b mod m$ ，其中 $a, b \leq m \leq 1 0^{18}$ 。

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

a \cdot b = ⎩ ⎨ ⎧ 0 2 \cdot \frac{a}{2} \cdot b 2 \cdot \frac{a - 1}{2} \cdot b + b if a = 0 if a > 0 and a even if a > 0 and a odd

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

高精度快速幂

前置技能：大整数乘法

洛谷 P1045 [NOIP 2003 普及组] 麦森数

给定整数 $P$ （ $1000 < P < 3100000$ ），计算 $2^{P} - 1$ 的位数与最后 $500$ 位数字（用十进制数表示），不足 $500$ 位时高位补 0。

代码实现

#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';
    }
  }
}

底数固定的预处理快速幂

当底数 $a$ 固定时，可以利用分块思想，用一定的时间预处理后用 $O (1)$ 的时间回答一次幂询问。这一算法也常称为光速幂。过程如下：

选定一个数 $s$ ，预处理出 $a^{0}, a^{1}, \dots, a^{s - 1}$ 与 $a^{0}, a^{s}, \dots, a^{⌊ p / s ⌋ s}$ 的值并存在两个数组里；
对于每一次询问 $a^{b}$ ，将 $b$ 拆分成 $⌊ b / s ⌋ s + (b mod s)$ ，则 $a^{b} = a^{⌊ b / s ⌋ s} \cdot a^{b mod s}$ ，就可以 $O (1)$ 求出答案。

假设指数 $b$ 的范围是 $[0, n]$ ，那么块长 $s$ 经常选择为 $n$ 或者与之相近的 $2$ 的幂次。选择 $n$ 可以获得最优的预处理复杂度 $O (n)$ ，而选择 $2$ 的幂次可以使用位运算简化计算。

特别地，对于模意义下幂的计算，底数 $a$ 相同隐含着模数 $m$ 也要相同这一要求。由于扩展欧拉定理，对于任意模数 $m$ ，预处理的指数范围上界为 $n = 2 φ (m)$ ；对于素模数 $p$ ，预处理的范围上界为 $n = p - 1$ 。这两种情形预处理的复杂度都是 $O (m)$ 。

参考代码

int a, mod, pow1[65536], pow2[65536];
 
void preproc() {
  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;
}

习题

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

Sean's Blog

Explorer

06.快速幂

引入

过程

迭代版本

递归版本

应用

模意义下取幂

计算斐波那契数

多次置换

加速几何中对点集的操作

定长路径计数

模意义下的整数乘法

高精度快速幂

底数固定的预处理快速幂

习题

目录

Graph View

反向链接