跳转至

快速冪

定義

快速冪，二進制取冪（Binary Exponentiation，也稱平方法），是一個在 \(\Theta(\log n)\) 的時間內計算 \(a^n\) 的小技巧，而暴力的計算需要 \(\Theta(n)\) 的時間。

這個技巧也常常用在非計算的場景，因為它可以應用在任何具有結合律的運算中。其中顯然的是它可以應用於模意義下取冪、矩陣冪等運算，我們接下來會討論。

解釋

計算 \(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\)。二進制取冪的想法是，我們將取冪的任務按照指數的 二進制表示 來分割成更小的任務。

過程

迭代版本

首先我們將 \(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{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} \]

因此為了計算 \(3^{13}\)，我們只需要將對應二進制位為 1 的整係數冪乘起來就行了：

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

將上述過程説得形式化一些，如果把 \(n\) 寫作二進制為 \((n_tn_{t-1}\cdots n_1n_0)_2\)，那麼有：

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

其中 \(n_i\in\{0,1\}\)。那麼就有

\[ \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} \]

根據上式我們發現，原問題被我們轉化成了形式相同的子問題的乘積，並且我們可以在常數時間內從 \(2^i\) 項推出 \(2^{i+1}\) 項。

這個算法的複雜度是 \(\Theta(\log n)\) 的，我們計算了 \(\Theta(\log n)\) 個 \(2^k\) 次冪的數，然後花費 \(\Theta(\log n)\) 的時間選擇二進制為 1 對應的冪來相乘。

遞歸版本

上述迭代版本中，由於 \(2^{i+1}\) 項依賴於 \(2^i\)，使得其轉換為遞歸版本比較困難（一方面需要返回一個額外的 \(a^{2^i}\)，對函數來説無法實現一個只返回計算結果的接口；另一方面則是必須從低位往高位計算，即從高位往低位調用，這也造成了遞歸實現的困擾），下面則提供遞歸版本的思路。

給定形式 \(n_{t\dots i} = (n_tn_{t-1}\cdots n_i)_2\)，即 \(n_{t\dots i}\) 表示將 \(n\) 的前 \(t - i + 1\) 位二進制位當作一個二進制數，則有如下變換：

\[ \begin{aligned} n &= n_{t\dots 0} \\ &= 2\times n_{t\dots 1} + n_0\\ &= 2\times (2\times n_{t\dots 2} + n_1) + n_0 \\ &\cdots \end{aligned} \]

那麼有：

\[ \begin{aligned} a^n &= a^{n_{t\dots 0}} \\ &= a^{2\times n_{t\dots 1} + n_0} = \left(a^{n_{t\dots 1}}\right)^2 a^{n_0} \\ &= \left(a^{2\times n_{t\dots 2} + n_1}\right)^2 a^{n_0} = \left(\left(a^{n_{t\dots 2}}\right)^2 a^{n_1}\right)^2 a^{n_0} \\ &\cdots \end{aligned} \]

如上所述，在遞歸時，對於不同的遞歸深度是相同的處理：\(a^{n_{t\dots i}} = (a^{n_{t\dots (i+1)}})^2a^{n_i}\)，即將當前遞歸的二進制數拆成兩部分：最低位在遞歸出來時乘上去，其餘部分則變成新的二進制數遞歸進入更深一層作相同的處理。

可以觀察到，每遞歸深入一層則二進制位減少一位，所以該算法的時間複雜度也為 \(\Theta(\log n)\)。

實現

首先我們可以直接按照上述遞歸方法實現：

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
2
3
4
5
6
7
8
def binpow(a, b):
    if b == 0:
        return 1
    res = binpow(a, b // 2)
    if (b % 2) == 1:
        return res * res * a
    else:
        return res * res

第二種實現方法是非遞歸式的。它在循環的過程中將二進制位為 1 時對應的冪累乘到答案中。儘管兩者的理論複雜度是相同的，但第二種在實踐過程中的速度是比第一種更快的，因為遞歸會花費一定的開銷。

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

1
2
3
4
5
6
7
8
def binpow(a, b):
    res = 1
    while b > 0:
        if (b & 1):
            res = res * a
        a = a * a
        b >>= 1
    return res

模板：Luogu P1226

應用

模意義下取冪

問題描述

計算 \(x^n\bmod m\)。

這是一個非常常見的應用，例如它可以用於計算模意義下的乘法逆元。

既然我們知道取模的運算不會干涉乘法運算，因此我們只需要在計算的過程中取模即可。

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

1
2
3
4
5
6
7
8
9
def binpow(a, b, m):
    a = a % m
    res = 1
    while b > 0:
        if (b & 1):
            res = res * a % m
        a = a * a % m
        b >>= 1
    return res

注意：根據費馬小定理，如果 \(m\) 是一個質數，我們可以計算 \(x^{n\bmod (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 \textit{length})\) 的時間內將這些變換應用到這個 \(n\) 個點，其中 \(\textit{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 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\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} \]

快速乘

但是 \(O(\log_2 m)\) 的「龜速乘」還是太慢了，這在很多對常數要求比較高的算法比如 Miller_Rabin 和 Pollard-Rho 中，就顯得不夠用了。所以我們要介紹一種可以處理模數在 long long 範圍內、不需要使用黑科技 __int128 的、複雜度為 \(O(1)\) 的「快速乘」。

我們發現：

\[ a\times b\bmod m=a\times b-\left\lfloor \dfrac{ab}m \right\rfloor\times m \]

我們巧妙運用 unsigned long long 的自然溢出：

\[ a\times b\bmod m=a\times b-\left\lfloor \dfrac{ab}m \right\rfloor\times m=\left(a\times b-\left\lfloor \dfrac{ab}m \right\rfloor\times m\right)\bmod 2^{64} \]

於是在算出 \(\left\lfloor\dfrac{ab}m\right\rfloor\) 後，兩邊的乘法和中間的減法部分都可以使用 unsigned long long 直接計算，現在我們只需要解決如何計算 \(\left\lfloor\dfrac {ab}m\right\rfloor\)。

我們考慮先使用 long double 算出 \(\dfrac am\) 再乘上 \(b\)。

既然使用了 long double，就無疑會有精度誤差。極端情況就是第一個有效數字（二進制下）在小數點後一位。在 x86-64 機器下，long double 將被解釋成 \(80\) 位拓展小數（即符號為 \(1\) 位，指數為 \(15\) 位，尾數為 \(64\) 位），所以 long double 最多能精確表示的有效位數為 \(64\)1。所以 \(\dfrac am\) 最差從第 \(65\) 位開始出錯，誤差範圍為 \(\left(-2^{-64},2^{64}\right)\)。乘上 \(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)\bmod m\)。

代碼實現如下：

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

高精度快速冪

前置技能

請先學習 高精度

例題【NOIP2003 普及組改編·麥森數】（原題在此）

題目大意：從文件中輸入 \(P\)（\(1000 < P < 3100000\)），計算 \(2^P−1\) 的最後 100 位數字（用十進制高精度數表示），不足 100 位時高位補 0。

代碼實現如下：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
#include <bits/stdc++.h>
using namespace std;
int a[505], b[505], t[505], i, j;

void 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);  //(2^(p/2))^2=2^p
  mult(b, b);  // 对b平方
  if (p % 2 == 1) mult(b, a);
}

int main() {
  int p;
  scanf("%d", &p);
  a[0] = 1;  // 记录a数组的位数
  a[1] = 2;  // 对2进行平方
  b[0] = 1;  // 记录b数组的位数
  b[1] = 1;  // 答案数组
  ksm(p);
  for (i = 100; i >= 1; i--) {
    if (i == 1) {
      printf("%d\n", b[i] - 1);  // 最后一位减1
    } else
      printf("%d", b[i]);
  }
}

同一底數與同一模數的預處理快速冪

在同一底數與同一模數的條件下，可以利用分塊思想，用一定的時間（一般是 \(O(\sqrt n)\)）預處理後用 \(O(1)\) 的時間回答一次冪詢問。

過程

  1. 選定一個數 \(s\)，預處理出 \(a^0\) 到 \(a^s\) 與 \(a^{0\cdot s}\) 到 \(a^{\lceil\frac ps\rceil\cdot s}\) 的值並存在一個（或兩個）數組裏；
  2. 對於每一次詢問 \(a^b\bmod p\)，將 \(b\) 拆分成 \(\left\lfloor\dfrac bs\right\rfloor\cdot s+b\bmod s\)，則 \(a^b=a^{\lfloor\frac bs\rfloor\cdot s}\times a^{b\bmod s}\)，可以 \(O(1)\) 求出答案。

關於這個數 \(s\) 的選擇，我們一般選擇 \(\sqrt p\) 或者一個大小適當的 \(2\) 的次冪（選擇 \(\sqrt p\) 可以使預處理較優，選擇 \(2\) 的次冪可以使用位運算優化/簡化計算）。

參考代碼 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
int pow1[65536], pow2[65536];

void preproc(int a, int mod) {
  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]; }

習題

本頁面部分內容譯自博文 Бинарное возведение в степень 與其英文翻譯版 Binary Exponentiation。其中俄文版版權協議為 Public Domain + Leave a Link；英文版版權協議為 CC-BY-SA 4.0。

參考資料與註釋

  1. 參見 C 語言小數表示法 - 維基百科 ↩

本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用