贰猹の小窝

【引自 OI Wiki】

因为以后可能会增加一点自己的理解，所以不打算每次都去查原帖，故在此直接做笔记和更改

---

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

计算 a 的 n 次方表示将 n 个 a 乘在一起：个 an=a×a⋯×a⏟n 个 a。然而当 a,n 太大的时侯，这种方法就不太适用了。不过我们知道：ab+c=ab⋅ac,a2b=ab⋅ab=(ab)2。二进制取幂的想法是，我们将取幂的任务按照指数的 二进制表示 来分割成更小的任务。

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

313=3(1101)2=38⋅34⋅31

因为 n 有 ⌊log2⁡n⌋+1 个二进制位，因此当我们知道了 a1,a2,a4,a8,…,a2⌊log2⁡n⌋ 后，我们只用计算 Θ(log⁡n) 次乘法就可以计算出 an。

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

(1)31=3(2)32=(31)2=32=9(3)34=(32)2=92=81(4)38=(34)2=812=6561

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

313=6561⋅81⋅3=1594323

将上述过程说得形式化一些，如果把 n 写作二进制为 (ntnt−1⋯n1n0)2，那么有：

n=nt2t+nt−12t−1+nt−22t−2+⋯+n121+n020

其中 ni∈0,1。那么就有

an=(ant2t+⋯+n020)=an020×an121×⋯×ant2t

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

这个算法的复杂度是 Θ(log⁡n) 的，我们计算了 Θ(log⁡n) 个 2k 次幂的数，然后花费 Θ(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

模意义下取幂

题目描述

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

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

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 是一个质数，我们可以计算 xnmod(m−1) 来加速算法过程。

计算斐波那契数

题目描述

根据斐波那契数列的递推式 Fn=Fn−1+Fn−2，我们可以构建一个 2×2 的矩阵来表示从 Fi,Fi+1 到 Fi+1,Fi+2 的变换。于是在计算这个矩阵的 n 次幂的时侯，我们使用快速幂的思想，可以在 Θ(log⁡n) 的时间内计算出结果。对于更多的细节参见 斐波那契数列。

题目描述

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

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

加速几何中对点集的操作

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

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

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

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

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

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

[a11a12a13a14a21a22a23a24a31a32a33a34a41a42a43a44]

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

[a11a12a13a14a21a22a23a24a31a32a33a34a41a42a43a44]⋅[xyz1]=[x′y′z′1]

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

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

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

   [1005010700190001]

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

   [10000050000500001]

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

   [10000cos⁡θsin⁡θ00−sin⁡θcos⁡θ00001]

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

定长路径计数

题目描述

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

模意义下大整数乘法

计算 a×bmodm,,,a,b≤m≤1018。

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

a⋅b={0if a=02⋅a2⋅bif a>0 and a even2⋅a−12⋅b+bif a>0 and a odd

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

我们发现：

a×bmodm=a×b−⌊abm⌋×m

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

a×bmodm=a×b−⌊abm⌋×m=(a×b−⌊abm⌋×m)mod264

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

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

既然使用了 long double，就无疑会有进度误差。极端情况就是第一个有效数字在小数点后一位。因为 sizeof(long double)=16，即 long double 的进度是 64 位有效数字。所以 ap 从第 65 位开始出错，误差范围为 (−2−64,264)。乘上 b 这个 64 位整数，误差范围为 (−0.5,0.5)，再加上 0.5 误差范围为 (0,1)，取整后误差范围位 0,1。于是乘上 −m 后，误差范围变成 0,−m，我们需要判断这两种情况。

因为 m 在 long long 范围内，所以如果计算结果 r 在 [0,m) 时，直接返回 r，否则返回 r+m，当然你也可以直接返回 (r+m)modm。

代码实现如下：

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

高精度快速幂

例题【NOIP2003 普及组改编・麦森数】（原题在此）

题目描述

代码实现如下：

#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]);
  }
}