快速幂and快速乘填空题

定义

快速幂，二进制取幂（Binary Exponentiation，也称平方法），是一个在  的时间内计算  的小技巧，而暴力的计算需要  的时间。

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

解释

计算  的  次方表示将  个  乘在一起：个。然而当  太大的时侯，这种方法就不太适用了。不过我们知道：。二进制取幂的想法是，我们将取幂的任务按照指数的 二进制表示 来分割成更小的任务。

过程

迭代版本

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

因为  有  个二进制位，因此当我们知道了  后，我们只用计算  次乘法就可以计算出 。

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

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

将上述过程说得形式化一些，如果把  写作二进制为 ，那么有：

其中 。那么就有

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

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

递归版本

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

给定形式 ，即  表示将  的前  位二进制位当作一个二进制数，则有如下变换：

那么有：

如上所述，在递归时，对于不同的递归深度是相同的处理：，即将当前递归的二进制数拆成两部分：最低位在递归出来时乘上去，其余部分则变成新的二进制数递归进入更深一层作相同的处理。

可以观察到，每递归深入一层则二进制位减少一位，所以该算法的时间复杂度也为 。

实现

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

C++Python

1 2 3 4 5 6 7 8

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

C++Python

1 2 3 4 5 6 7 8 9

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

应用

模意义下取幂

问题描述

计算 。

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

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

C++Python

1 2 3 4 5 6 7 8 9 10

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

注意：根据费马小定理，如果  是一个质数，我们可以计算  来加速算法过程。

计算斐波那契数

问题描述

计算斐波那契数列第  项 。

根据斐波那契数列的递推式 ，我们可以构建一个  的矩阵来表示从  到  的变换。于是在计算这个矩阵的  次幂的时侯，我们使用快速幂的思想，可以在  的时间内计算出结果。对于更多的细节参见 斐波那契数列，矩阵快速幂的实现参见 矩阵加速递推 中的实现。

多次置换

问题描述

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

简单地把这个置换取  次幂，然后把它应用到序列  上即可。时间复杂度是  的。

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

加速几何中对点集的操作

引入

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

1. 沿某个向量移动点的位置（Shift）。
2. 按比例缩放这个点的坐标（Scale）。
3. 绕某个坐标轴旋转（Rotate）。

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

解释

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

1. Shift 操作：将每一维的坐标分别加上一个常量；
2. Scale 操作：把每一维坐标分别乘上一个常量；
3. Rotate 操作：这个有点复杂，我们不打算深入探究，不过我们仍然可以使用一个线性组合来表示新的坐标。

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

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

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

过程

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

1. Shift 操作：让  坐标方向的位移为 ， 坐标的位移为 ， 坐标的位移为 ：

   

2. Scale 操作：把  坐标拉伸 10 倍， 坐标拉伸 5 倍：

   

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

   

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

定长路径计数

问题描述

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

我们把该图的邻接矩阵 M 取 k 次幂，那么  就表示从  到  长度为  的路径的数目。该算法的复杂度是 。有关该算法的细节请参见 矩阵 页面。

模意义下大整数乘法

计算 。

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

快速乘

但是  的「龟速乘」还是太慢了，这在很多对常数要求比较高的算法比如 Miller_Rabin 和 Pollard-Rho 中，就显得不够用了。所以我们要介绍一种可以处理模数在 long long 范围内、不需要使用黑科技 __int128 的、复杂度为  的「快速乘」。

我们发现：

我们巧妙运用 unsigned long long 的自然溢出：

于是在算出  后，两边的乘法和中间的减法部分都可以使用 unsigned long long 直接计算，现在我们只需要解决如何计算 。

我们考虑先使用 long double 算出  再乘上 。

既然使用了 long double，就无疑会有精度误差。极端情况就是第一个有效数字（二进制下）在小数点后一位。在 x86-64 机器下，long double 将被解释成  位拓展小数（即符号为  位，指数为  位，尾数为  位），所以 long double 最多能精确表示的有效位数为 1。所以  最差从第  位开始出错，误差范围为 。乘上  这个  位整数，误差范围为 ，再加上  误差范围为 ，取整后误差范围位 。于是乘上  后，误差范围变成 ，我们需要判断这两种情况。

因为  在 long long 范围内，所以如果计算结果  在  时，直接返回 ，否则返回 ，当然你也可以直接返回 。

代码实现如下：

1 2 3 4 5 6 7

long long binmul(long long a, long long b, long long m) { unsigned long long c = (unsigned long long)a * b - (unsigned long long)((long double)a / m * b + 0.5L) * m; if (c < m) return c; return c + m; }