Merge pull request #1416 from sshwy/quick-pow

快速幂，是一种求 $a^b \bmod p$ 的方法，得益于将指数按二进制拆开的思想。

快速幂，二进制取幂（Binary Exponentiation，也称平方法），是一个在 $\Theta(\log_2n)$ 的时间内计算 $a^n$ 的小技巧，而暴力的计算需要 $\Theta(n)$ 的时间。而这个技巧也常常用在非计算的场景，因为它可以应用在任何具有结合律的运算中。其中显然的是它可以应用于模意义下取幂、矩阵幂等运算，我们接下来会讨论。

事实上，根据模运算的性质， $a \times b \bmod p = ((a \bmod p) \times b) \bmod p$ 。那么我们也可以把 $a^b \mod p$ 分解成一系列比较小的数的乘积。

## 算法描述

如果把 $b$ 写作二进制为 $a_ta_{t-1} \cdots a_1a_0$ ，那么有：

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

$$

b = a_t2^t + a_{t-1}2^{t-1} + a_{t-2}2^{t-2} + \cdots + a_12^1 + a_02^0

$$

，其中 $a_i$ 是 0 或者 1。

那么就有

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

$$

\begin{aligned}

a^b \bmod p & = (a^{a_t 2^t + \cdots + a_0 2^0}) \bmod p \\\\

& = (..(a^{a_0 2^0} \bmod p) \times \cdots \times a^{a_52^5}) \bmod p

\end{aligned}

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}}$ 后，我们只用计算 $\Theta(\log_2n)$ 次乘法就可以计算出 $a^n$ 。

最重要的是，我们注意到， $a^{2^{i+1}} \bmod c = (a^{2^i})^2 \bmod c$ ，可以在常数时间内从 $2^i$ 项推出 $2^{i+1}$ 项。于是，原问题总的复杂度就是 $O(\log b)$ 

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

在算法竞赛中，快速幂的思想不仅用于整数乘法，也可用于大整数加法，矩阵幂运算等场合中。

如果你看不懂，那就简单点说吧。

$$

\begin{align}

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{align}

$$

举个栗子， $a^{10}$ 等价于下面的式子：

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

 $a \times a \times a \times a \times a \times a \times a \times a \times a \times a$ 

$$

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

$$

通过观察我们不难发现， $a^{10}$ 可以转化成 $(a \times a)^{5}$ 

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

 $\left(a \times a \right) \times\left(a \times a \right) \times \left(a \times a \right) \times \left(a \times a \right) \times \left(a \times a \right)$ 

$$

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

$$

这时，再进行分解，我们假设 $a' =a \times a$ ，原式就成了

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

$$

a'\times a'\times a'\times a'\times a'

\begin{aligned}

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_t2^t}

\end{aligned}

$$

可是我们发现，a 不能正好分完，于是我们单独拎出来一个 a'，就转化成了：

 $\left (a' \times a'\right) \times\left (a' \times a'\right) \times a'$ 

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

如此重复下去即可，终止条件：

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

 $a^0=1$ 和 $a^1=a$ 

## 代码实现

## 实现代码

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

注意，这种方法能实现的问题比较单调，不可以解决大整数加法，矩阵幂运算。

```cpp

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;

}

```

### 非递归版

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

```cpp

int quickPow(int a, int b, int c) {

  // calculates a^b mod c

  int res = 1, bas = a;

  while (b) {

    if (b & 1) res = (LL)res * bas % c;

    // Transform to long long in case of overflow.

    bas = bas * bas % c;

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;

}

```

### 递归版

??? note "例题"

    做一做 [Luogu P1226](https://www.luogu.org/problemnew/show/P1226)

## 应用

### 模意义下取幂

???+note "问题描述"

    计算 $x^n\bmod m$。

这是一个非常常见的应用，例如它可以用于计算模意义下的乘法逆元。

既然我们知道取模的运算不会干涉乘法运算，因此我们只需要在计算的过程中取模即可。

```cpp

long long qpow(long long a, long long b, long long p) {

  if (b == 0) return 1 % p;

  if (b == 1) return a % p;

  if (b % 2 == 0) {

    long long t = a * a % p;

    return qpow(t, b / 2, p);

  } else {

    long long t = a * a % p;

    return (qpow(t, b / 2, p) * a) % p;

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;

}

```

??? note "例题"

注意：根据费马小定理，如果 $m$ 是一个质数，我们可以计算 $x^{n\bmod (m-1)}$ 来加速算法过程。

    做一做[Luogu P1226](https://www.luogu.org/problemnew/show/P1226)

### 计算斐波那契数

???+note "问题描述"

    计算斐波那契数列第 $n$ 项 $F_n$。

根据斐波那契数列的递推式 $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)$ 的时间内计算出结果。对于更多的细节参见[OI-wiki 斐波那契数列](/math/fibonacci/)。

### 多次置换

???+note "问题描述"

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

简单地把这个置换取 $k$ 次幂，然后把它应用到序列 $n$ 上即可。时间复杂度是 $O(n \log k)$ 的。

注意：给这个置换建图，然后在每一个环上分别做 $k$ 次幂（事实上做一下 $k$ 对环长取模的运算即可）可以取得更高效的算法，达到 $O(n)$ 的复杂度。

### 加速几何中对点集的操作

## 高精度快速幂

> 三维空间中， $n$ 个点 $p_i$ ，要求将 $m$ 个操作都应用于这些点。包含 3 种操作：

>

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

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

> 3.  绕某个坐标轴旋转（Rotate）。

>

> 还有一个特殊的操作，就是将一个操作序列重复 $k$ 次（Loop），这个序列中也可能有 Loop 操作（Loop 操作可以嵌套）。现在要求你在低于 $O(n \cdot length)$ 的时间内将这些变换应用到这个 $n$ 个点，其中 $length$ 表示把所有的 Loop 操作展开后的操作序列的长度。

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

1.  Shift 操作：将每一维的坐标分别加上一个常量；

2.  Scale 操作：把每一维坐标分别乘上一个常量；

3.  Rotate 操作：这个有点复杂，我们不打算深入探究，不过我们仍然可以使用一个线性组合来表示新的坐标。

可以看到，每一个变换可以被表示为对坐标的线性运算，因此，一个变换可以用一个 $4\times 4$ 的矩阵来表示：

$$

\begin{bmatrix}

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{bmatrix}

$$

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

$$

\begin{bmatrix}

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{bmatrix}\cdot

\begin{bmatrix} x \\ y \\ z \\ 1 \end{bmatrix}

 = \begin{bmatrix} x' \\ y' \\ z' \\ 1 \end{bmatrix}

$$

你可能会问，为什么一个三维坐标会多一个 1 出来？原因在于，如果没有这个多出来的 1，我们没法使用矩阵的线性变换来描述 Shift 操作。

接下来举一些简单的例子来说明我们的思路：

1.  Shift 操作：让 $x$ 坐标方向的位移为 $5$ ， $y$ 坐标的位移为 $7$ ， $z$ 坐标的位移为 $9$ ：

    $$

    \begin{bmatrix}

    1 & 0 & 0 & 5 \\

    0 & 1 & 0 & 7 \\

    0 & 0 & 1 & 9 \\

    0 & 0 & 0 & 1 \\

    \end{bmatrix}

    $$

2.  Scale 操作：把 $x$ 坐标拉伸 10 倍， $y,z$ 坐标拉伸 5 倍：

    $$

    \begin{bmatrix}

    10 & 0 & 0 & 0 \\

    0 & 5 & 0 & 0 \\

    0 & 0 & 5 & 0 \\

    0 & 0 & 0 & 1 \\

    \end{bmatrix}

    $$

3.  Rotate 操作：绕 $x$ 轴旋转 $\theta$ 弧度，遵循右手定则（逆时针方向）

    $$

    \begin{bmatrix}

    1 & 0 & 0 & 0 \\

    0 & \cos \theta & \sin \theta & 0 \\

    0 & -\sin \theta & \cos \theta & 0 \\

    0 & 0 & 0 & 1 \\

    \end{bmatrix}

    $$

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

### 定长路径计数

???+note "问题描述"

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

我们把该图的邻接矩阵 M 取 k 次幂，那么 $M_{i,j}$ 就表示从 $i$ 到 $j$ 长度为 $k$ 的路径的数目。该算法的复杂度是 $O(n^3 \log_2 k)$ 。有关该算法的细节参见[Number of paths of fixed length/Shortest paths of fixed length](https://cp-algorithms.com/graph/fixed_length_paths.html)。

### 模意义下大整数乘法

> 计算 $a\times b\bmod m,\,\,a,b\le m\le 10^{18}$ 。

与二进制取幂的思想一样，这次我们将其中的一个乘数表示为若干个 2 的整数次幂的和的形式。因为在对一个数做乘 2 并取模的运算的时侯，我们可以转化为加减操作防止溢出。这样仍可以在 $O (\log_2 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}

$$

注意：你也可以利用双精度浮点数在常数时间内计算大整数乘法。因为 $a\times b\bmod m=a\times b-\left\lfloor\frac{a\times b}{m}\right\rfloor m$ 。由于 $a,b<m$ ，因此 $\left\lfloor\frac{a\times b}{m}\right\rfloor<m$ ，于是可以用双精度浮点数计算这个分式。作差的时侯直接自然溢出。因为两者的差是一定小于 $m$ 的，我们只关心低位。这样再调整一下正负性就行了。更多信息参见[这里](https://cs.stackexchange.com/questions/77016/modular-multiplication)。

### 高精度快速幂

??? note "前置技能"

     请先学习[高精度](https://oi-wiki.org/math/bignum/)

??? note " 例题【NOIP2003 普及组改编·麦森数】（[原题在此](https://www.luogu.org/problemnew/show/P1045)）"

???+note " 例题【NOIP2003 普及组改编·麦森数】（[原题在此](https://www.luogu.org/problemnew/show/P1045)）"

     题目大意：从文件中输入 P（1000<P<3100000），计算 $2^P−1$ 的最后 100 位数字（用十进制高精度数表示），不足 100 位时高位补 0。

### 例题代码

代码实现如下：

```cpp

#include <bits/stdc++.h>

  }

}

```

## 习题

-   [UVa 1230 - MODEX](http://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&category=24&page=show_problem&problem=3671)

-   [UVa 374 - Big Mod](http://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&category=24&page=show_problem&problem=310)

-   [UVa 11029 - Leading and Trailing](https://uva.onlinejudge.org/index.php?option=onlinejudge&page=show_problem&problem=1970)

-   [Codeforces - Parking Lot](http://codeforces.com/problemset/problem/630/I)

-   [SPOJ - The last digit](http://www.spoj.com/problems/LASTDIG/)

-   [SPOJ - Locker](http://www.spoj.com/problems/LOCKER/)

-   [LA - 3722 Jewel-eating Monsters](https://icpcarchive.ecs.baylor.edu/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=1723)

-   [SPOJ - Just add it](http://www.spoj.com/problems/ZSUM/)

     **本页面部分内容译自博文[Бинарное возведение в степень](http://e-maxx.ru/algo/binary_pow)与其英文翻译版[Binary Exponentiation](https://cp-algorithms.com/algebra/binary-exp.html)。其中俄文版版权协议为 Public Domain + Leave a Link；英文版版权协议为 CC-BY-SA 4.0。**