【Math】快速幂

百科名片

快速幂，二进制取幂（Binary Exponentiation，也称平方法），是一个在 $\Theta (\log n)$ 的时间内计算 $a^n$ 的小技巧，而暴力的计算需要 $\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 }}$ 后，我们只用计算 $\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}$ 项。

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

代码实现

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

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

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

模板： Luogu-P1226

应用

模意义下取幂

问题描述

计算 $x^n\mathrm{m}\mathrm{o}\mathrm{d}m$ 。

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

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

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$ 是一个质数，我们可以计算 $x^{n\mathrm{m}\mathrm{o}\mathrm{d}(m-1)}$ 来加速算法过程。

计算斐波那契数

问题描述

计算斐波那契数列第 $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)$ 的时间内计算出结果。

多次置换

问题描述

给你一个长度为 $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$ 的矩阵来表示：

$$\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 出来？原因在于，如果没有这个多出来的 1，我们没法使用矩阵的线性变换来描述 Shift 操作。

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

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

   $$\left[\begin{array}{cccc}1& 0& 0& 5\\ 0& 1& 0& 7\\ 0& 0& 1& 9\\ 0& 0& 0& 1\end{array}\right]$$

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

   $$\left[\begin{array}{cccc}10& 0& 0& 0\\ 0& 5& 0& 0\\ 0& 0& 5& 0\\ 0& 0& 0& 1\end{array}\right]$$

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

   $$\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& \cos \theta & \sin \theta & 0\\ 0& -\sin \theta & \cos \theta & 0\\ 0& 0& 0& 1\end{array}\right]$$

现在，每一种操作都被表示为了一个矩阵，变换序列可以用矩阵的乘积来表示，而一个 Loop 操作相当于取一个矩阵的 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\times b\mathrm{m}\mathrm{o}\mathrm{d}m,a,b\le m\le 10^{18}$ 。

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

注意：你也可以利用双精度浮点数在常数时间内计算大整数乘法。因为 $a\times b\mathrm{m}\mathrm{o}\mathrm{d}m=a\times b-\lfloor \frac{a\times b}{m}\rfloor m$ 。由于 $a,b<m$ ，因此 $\lfloor \frac{a\times b}{m}\rfloor <m$ ，于是可以用双精度浮点数计算这个分式。作差的时侯直接自然溢出。因为两者的差是一定小于 $m$ 的，我们只关心低位。这样再调整一下正负性就行了。

高精度快速幂

前置技能

请先学习 高精度

例题

【NOIP2003 普及组改编·麦森数】（原题在此）
题目大意：从文件中输入 P（1000<P<3100000），计算 $2^P-1$ 的最后 100 位数字（用十进制高精度数表示），不足 100 位时高位补 0。

代码实现如下：

#include <bits/stdc++.h>
using namespace std;
int a[505], b[505], t[505], i, j;

int mult(int x[], int y[]) { // 高精度乘法
	memset(t, 0, sizeof(t));
	for (i = 1; i <= x[0]; i++) {
		for (j = 1; j <= y[0]; j++) {
			if (i + j - 1 > 100) 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 ksm(int p) { // 快速幂
	if (p == 1) {
		memcpy(b, a, sizeof(b));
		return;
	}
	ksm(p / 2);
	mult(b, b);
	if (p % 2 == 1) mult(b, a);
}

int main() {
	int p;
	scanf("%d", &p);
	a[0] = 1;
	a[1] = 2;
	b[0] = 1;
	b[1] = 1;
	ksm(p);
	for (i = 100; i >= 1; i--) {
		if (i == 1) {
			printf("%d\n", b[i] - 1);
		} else printf("%d", b[i]);
	}
	return 0;
}

相关练习

• UVa-1230 MODEX
  原题 vjudge

• UVa-374 Big Mod
  原题 vjudge

• UVa-11029 Leading and Trailing
  原题 vjudge

• CodeForces-630I Parking Lot
  原题 vjudge

• SPOJ The last digit

• 原题 vjudge

• SPOJ Magic of the locker
  原题 vjudge

• LA-3722 Jewel-eating Monsters
  原题 vjudge

• SPOJ Just Add It
  原题 vjudge