Lũy thừa nhị phân

引入

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

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

过程

计算 \(a\) 的 \(n\) 次方表示将 \(n\) 个 \(a\) 乘在一起：\(a^{n} = \underbrace{a \times a \cdots \times a}_{n\text{ 个 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}\) 的计算过程如下：

\[ \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, \\ 3^{13} &= 6561 \times 81 \times 3 = 1594323. \end{aligned} \]

过程中，只进行了 \(5\) 次乘法运算．

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

迭代版本

设 \(n\) 的二进制表示为 \((n_tn_{t-1}\cdots n_1n_0)_2\)，也就是说，有

\[ n = n_t2^t + n_{t-1}2^{t-1} + \cdots + n_12^1 + n_02^0, \]

其中，\(n_i\in\{0,1\}\)．那么，就有

\[ \begin{aligned} a^n & = a^{n_t2^t + n_{t-1}2^{t-1} + \cdots + n_12^1 + n_02^0}\\ & = a^{n_0 2^0} \times a^{n_1 2^1}\times \cdots \times a^{n_{t-1}2^{t-1}} \times a^{n_t2^t}. \end{aligned} \]

注意，只有 \(n_i=1\) 的项才会真正出现在乘积的计算中．

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

伪代码如下：

\[ \begin{array}{l} \textbf{Algorithm }\text{FastPow}(a, n): \\ \textbf{Input. }\text{Base }a\text{ and exponent }n.\\ \textbf{Output. }\text{Power }a^n.\\ \textbf{Method.}\\ \begin{array}{ll} 1 & \textit{result}\gets\mathrm{Id}\\ 2 & \textbf{while }n > 0\textbf{ do}\\ 3 & \qquad \textbf{if }n \bmod 2 = 1\textbf{ then}\\ 4 & \qquad \qquad \textit{result} \gets \textit{result}\cdot a\\ 5 & \qquad \textbf{end if}\\ 6 & \qquad a \gets a \cdot a\\ 7 & \qquad n \gets n / 2\\ 8 & \textbf{end while}\\ 9 & \textbf{return }\textit{result} \end{array} \end{array} \]

利用这一方法计算快速幂，需要进行 \(\Theta(\log n)\) 次乘法运算．

递归版本

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

\[ (n_tn_{t-1}\cdots n_1n_0)_2 = 2 \times (n_tn_{t-1}\cdots n_1)_2 + n_0. \]

因此，幂次 \(a^n\) 可以递归地计算为

\[ a^n = \begin{cases} 1, & n = 0,\\ (a^{\lfloor n/2\rfloor})^2, & n > 0 \text{ and }n\text{ is even},\\ (a^{\lfloor n/2\rfloor})^2\cdot a, & n > 0 \text{ and }n\text{ is odd}.\\ \end{cases} \]

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

伪代码如下：

\[ \begin{array}{l} \textbf{Algorithm }\text{FastPow}(a, n): \\ \textbf{Input. }\text{Base }a\text{ and exponent }n.\\ \textbf{Output. }\text{Power }a^n.\\ \textbf{Method.}\\ \begin{array}{ll} 1 & \textbf{if }n = 0\textbf{ then}\\ 2 & \qquad \textbf{return }\mathrm{Id}\\ 3 & \textbf{end if}\\ 4 & \textit{result} \gets \text{FastPow}(a, n / 2) \\ 5 & \textbf{if }n\bmod 2 = 0\textbf{ then}\\ 6 & \qquad \textbf{return }\textit{result}\cdot\textit{result}\\ 7 & \textbf{else}\\ 8 & \qquad \textbf{return }\textit{result}\cdot\textit{result}\cdot a\\ 9 & \textbf{end if} \end{array} \end{array} \]

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

应用

模意义下取幂

洛谷 P1226【模板】快速幂

给定三个整数 \(a,b,p\)，求 \(a^b\bmod p\)．其中 \(p\ge 2\)．

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

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

参考实现

C++Python

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

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

参考实现

C++Python

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

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

注意

模数通常情况下大于 \(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\) 次幂的时侯，我们使用快速幂的思想，可以在 \(\Theta(\log n)\) 的时间内计算出结果．对于更多的细节参见斐波那契数列，矩阵快速幂的实现参见矩阵加速递推中的实现．

多次置换

问题描述

给你一个长度为 \(n\) 的序列和一个置换，把这个序列置换 \(k\) 次．

简单地把这个置换取 \(k\) 次幂，然后把它应用到序列上即可．时间复杂度为 \(O(n \log 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 \log k)\) 的时间计算出整个变换序列最终形成的矩阵．最后将它应用到 \(n\) 个点上，总复杂度 \(O(n + m \log k)\)．

定长路径计数

问题描述

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

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

模意义下的整数乘法

问题描述

给定非负整数 \(a,b\) 和正整数 \(m\)，计算 \(a\times b\bmod m\)，其中 \(a,b\le m\le 10^{18}\)．

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

\[ a \cdot b = \begin{cases} 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{cases} \]

但在实际使用中，此方法由于引入了更大的计算复杂度导致时间效率不优．实际编程中通常利用快速乘来进行模数范围在 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,\cdots,a^{s-1}\) 与 \(a^0,a^s,\cdots,a^{\lfloor p/s\rfloor s}\) 的值并存在两个数组里；
对于每一次询问 \(a^b\)，将 \(b\) 拆分成 \(\lfloor b/s\rfloor s+(b\bmod s)\)，则 \(a^b=a^{\lfloor b/s\rfloor s}\cdot a^{b\bmod s}\)，就可以 \(O(1)\) 求出答案．

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

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

Last updated on this page:, Update history
Found an error? Want to help improve? Edit this page on GitHub!
Contributors to this page:OI-wiki
All content on this page is provided under the terms of the CC BY-SA 4.0 and SATA license, additional terms may apply