Bỏ qua

Lũy thừa nhị phân

Giới thiệu

Lũy thừa nhanh (fast exponentiation), còn gọi là nhị phân lũy thừa (binary exponentiation) hay bình phương lũy thừa (exponentiation by squaring), là một kỹ thuật tính \(a^n\) trong \(\Theta(\log n)\) thời gian, trong khi cách vét cạn cần \(\Theta(n)\).

Kỹ thuật này áp dụng cho mọi bối cảnh mà phép nhân của \(a\) thỏa mãn tính kết hợp, như lũy thừa modulo, lũy thừa ma trận, v.v., xem thêm ở mục Ứng dụng.

Quy trình

Tính \(a^n\) tức là nhân \(n\) lần \(a\): \(a^{n} = \underbrace{a \times a \cdots \times a}_{n\text{ lần a}}\). Tuy nhiên khi \(n\) lớn hoặc phép nhân tốn kém, cách này không phù hợp. Ý tưởng của nhị phân lũy thừa là chia nhỏ bài toán theo biểu diễn nhị phân của số mũ.

Ví dụ

Cần tính \(3^{13}\). Nếu nhân liên tiếp thì cần \(13-1=12\) phép nhân. Nhưng vì

\[ 3^{13} = 3^{(1101)_2} = 3^8 \times 3^4 \times 3^1, \]

nên chỉ cần tính nhanh \(3^{1},3^{2},3^{4},3^{8}\) rồi nhân lại, mất \(2\) phép nhân để ra \(3^{13}\). Ta chỉ cần cách tính nhanh dãy các lũy thừa \(3^{2^k}\). Điều này dễ vì mỗi phần tử (trừ phần tử đầu) là bình phương của phần tử trước.

Quá trình tính \(3^{13}\):

\[ \begin{aligned} 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, \\ 3^{13} &= 6561 \times 81 \times 3 = 1594323. \end{aligned} \]

Quá trình chỉ thực hiện \(5\) phép nhân.

Đó là ý tưởng cơ bản. Có hai cách cài đặt thường gặp.

Phiên bản lặp

Gọi biểu diễn nhị phân của \(n\) là \((n_tn_{t-1}\cdots n_1n_0)_2\), tức

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

với \(n_i\in\{0,1\}\). Khi đó

\[ \begin{aligned} a^n & = a^{n_t2^t + n_{t-1}2^{t-1} + \cdots + n_12^1 + n_02^0}\\ & = a^{n_0 2^0} \times a^{n_1 2^1}\times \cdots \times a^{n_{t-1}2^{t-1}} \times a^{n_t2^t}. \end{aligned} \]

Chỉ các hạng với \(n_i=1\) mới xuất hiện trong tích.

Từ biểu thức này, ta có thể tính \(a^{2^k}\) trong \(\Theta(\log n)\) thời gian, rồi nhân các lũy thừa tương ứng với bit bằng \(1\) để ra kết quả trong \(\Theta(\log n)\) thời gian. Đó là phiên bản lặp.

Mã giả:

\[ \begin{array}{l} \textbf{Algorithm }\text{FastPow}(a, n): \\ \textbf{Input. }\text{Base }a\text{ and exponent }n.\\ \textbf{Output. }\text{Power }a^n.\\ \textbf{Method.}\\ \begin{array}{ll} 1 & \textit{result}\gets\mathrm{Id}\\ 2 & \textbf{while }n > 0\textbf{ do}\\ 3 & \qquad \textbf{if }n \bmod 2 = 1\textbf{ then}\\ 4 & \qquad \qquad \textit{result} \gets \textit{result}\cdot a\\ 5 & \qquad \textbf{end if}\\ 6 & \qquad a \gets a \cdot a\\ 7 & \qquad n \gets n / 2\\ 8 & \textbf{end while}\\ 9 & \textbf{return }\textit{result} \end{array} \end{array} \]

Cách này cần \(\Theta(\log n)\) phép nhân.

Phiên bản đệ quy

Quy trình cũng có thể viết đệ quy. Lưu ý biểu diễn nhị phân của \(n\) có thể viết:

\[ (n_tn_{t-1}\cdots n_1n_0)_2 = 2 \times (n_tn_{t-1}\cdots n_1)_2 + n_0. \]

Do đó

\[ a^n = \begin{cases} 1, & n = 0,\\ (a^{\lfloor n/2\rfloor})^2, & n > 0 \text{ và }n\text{ chẵn},\\ (a^{\lfloor n/2\rfloor})^2\cdot a, & n > 0 \text{ và }n\text{ lẻ}.\\ \end{cases} \]

Mã giả:

\[ \begin{array}{l} \textbf{Algorithm }\text{FastPow}(a, n): \\ \textbf{Input. }\text{Base }a\text{ and exponent }n.\\ \textbf{Output. }\text{Power }a^n.\\ \textbf{Method.}\\ \begin{array}{ll} 1 & \textbf{if }n = 0\textbf{ then}\\ 2 & \qquad \textbf{return }\mathrm{Id}\\ 3 & \textbf{end if}\\ 4 & \textit{result} \gets \text{FastPow}(a, n / 2) \\ 5 & \textbf{if }n\bmod 2 = 0\textbf{ then}\\ 6 & \qquad \textbf{return }\textit{result}\cdot\textit{result}\\ 7 & \textbf{else}\\ 8 & \qquad \textbf{return }\textit{result}\cdot\textit{result}\cdot a\\ 9 & \textbf{end if} \end{array} \end{array} \]

Cách này đệ quy \(\Theta(\log n)\) lần và cũng cần \(\Theta(\log n)\) phép nhân. Do overhead đệ quy, thực tế phiên bản lặp thường nhanh hơn.

Ứng dụng

Lũy thừa modulo

洛谷 P1226【模板】快速幂

Cho ba số nguyên \(a,b,p\), tính \(a^b\bmod p\) với \(p\ge 2\).

Đây là ứng dụng rất phổ biến, ví dụ để tính nghịch đảo nhân modulo. Vì phép lấy modulo không ảnh hưởng đến phép nhân, chỉ cần lấy modulo trong quá trình tính.

Ta có thể cài đặt trực tiếp theo bản đệ quy:

Cài đặt tham khảo 
1
2
3
4
5
6
7
8
long long binpow(long long a, long long b, long long p) {
  if (b == 0) return 1;
  long long res = binpow(a, b / 2, p);
  if (b % 2)
    return res * res % p * a % p;
  else
    return res * res % p;
}

1
2
3
4
5
6
7
8
def binpow(a, b, p):
    if b == 0:
        return 1
    res = binpow(a, b // 2, p)
    if (b % 2) == 1:
        return res * res * a % p
    else:
        return res * res % p

Cách thứ hai là không đệ quy. Trong vòng lặp, ta nhân vào kết quả khi bit nhị phân bằng 1. Dù độ phức tạp lý thuyết giống nhau, cách này thường nhanh hơn vì không có overhead đệ quy.

Cài đặt tham khảo 
1
2
3
4
5
6
7
8
9
long long binpow(long long a, long long b, long long p) {
  long long res = 1;
  while (b > 0) {
    if (b & 1) res = res * a % p;
    a = a * a % p;
    b >>= 1;
  }
  return res;
}

1
2
3
4
5
6
7
8
def binpow(a, b, p):
    res = 1
    while b > 0:
        if b & 1:
            res = res * a % p
        a = a * a % p
        b >>= 1
    return res
Lưu ý
  • Thông thường mô-đun lớn hơn \(1\). Trong trường hợp rất đặc biệt \(p=1\), cần xét riêng trường hợp \(b=0\).
  • Khi số mũ rất lớn, cần dùng định lý Euler mở rộng để giảm mũ.

Tính số Fibonacci

Theo truy hồi \(F_n = F_{n-1} + F_{n-2}\), ta có thể xây dựng ma trận \(2\times 2\) biểu diễn biến đổi từ \(F_i,F_{i+1}\) sang \(F_{i+1},F_{i+2}\). Khi tính ma trận này mũ \(n\), áp dụng lũy thừa nhanh sẽ cho kết quả trong \(\Theta(\log n)\). Xem thêm dãy Fibonacci và tăng tốc truy hồi bằng ma trận.

Lặp hoán vị nhiều lần

Mô tả bài toán

Cho một dãy độ dài \(n\) và một hoán vị, áp dụng hoán vị đó \(k\) lần.

Chỉ cần lấy lũy thừa \(k\) của hoán vị rồi áp dụng lên dãy. Độ phức tạp \(O(n \log k)\). Xem thêm hợp thành hoán vị.

Lưu ý

Xây dựng đồ thị của hoán vị và xử lý từng chu trình với \(k\) lần (tương đương \(k\) modulo độ dài chu trình) sẽ giải trong \(O(n)\).

Tăng tốc thao tác hình học trên tập điểm

HDU 4087 A Letter to Programmers

Cho \(n\) điểm \(p_i\) trong không gian 3D, áp dụng \(m\) thao tác, gồm:

  1. Dời điểm theo một vector (Shift).
  2. Co/giãn theo tỉ lệ (Scale).
  3. Quay quanh một đường thẳng (Rotate).

Còn có thao tác đặc biệt lặp lại một chuỗi thao tác \(k\) lần (Repeat), Repeat có thể lồng nhau. Yêu cầu tọa độ cuối cùng.

Tham khảo Vector và ma trận. Mỗi thao tác là một ma trận biến đổi, chuỗi thao tác là tích ma trận. Repeat tương đương lấy lũy thừa ma trận. Nhờ đó tính ma trận tổng hợp trong \(O(m \log k)\), rồi áp dụng lên \(n\) điểm, tổng \(O(n + m \log k)\).

Đếm đường đi độ dài cố định

Mô tả bài toán

Cho đồ thị có hướng (trọng số cạnh bằng 1), đếm số đường đi độ dài \(k\) từ \(u\) đến \(v\) cho mọi cặp.

Lấy ma trận kề \(M\) mũ \(k\), khi đó \(M_{i,j}\) là số đường đi độ dài \(k\) từ \(i\) đến \(j\). Độ phức tạp \(O(n^3 \log k)\). Xem Ma trận.

Nhân số nguyên modulo

Mô tả bài toán

Cho \(a,b\) không âm và \(m\) dương, tính \(a\times b\bmod m\), với \(a,b\le m\le 10^{18}\).

Tương tự lũy thừa nhanh, ta biểu diễn một thừa số thành tổng các lũy thừa của 2. Vì nhân 2 và lấy modulo có thể chuyển thành cộng/trừ để tránh tràn, nên bài toán giải được trong \(O(\log m)\). Công thức đệ quy:

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

Tuy nhiên trong thực tế, cách này do phức tạp lớn hơn nên kém hiệu quả. Thường dùng nhân nhanh khi mô-đun trong long long.

Lũy thừa nhanh với số lớn

Kỹ năng trước: nhân số lớn

洛谷 P1045 [NOIP 2003 普及组] 麦森数

Cho \(P\) (\(1000 < P < 3100000\)), tính số chữ số của \(2^P−1\) và 500 chữ số cuối (thập phân), nếu thiếu thì thêm 0 ở đầu.

Cài đặt 
 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
46
47
48
49
50
51
52
53
#include <cmath>
#include <cstring>
#include <iostream>
using namespace std;

const int M = 500;

int a[505], b[505], t[505];

// 大整数乘法
void mult(int x[], int y[]) {
  memset(t, 0, sizeof(t));
  for (int i = 1; i <= x[0]; i++) {
    for (int j = 1; j <= y[0]; j++) {
      if (i + j - 1 > M) 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 binpow(int p) {
  if (p == 1) {
    memcpy(b, a, sizeof(b));
    return;
  }
  binpow(p / 2);  // (2^(p/2))^2=2^p
  mult(b, b);     // 对 b 平方
  if (p % 2 == 1) mult(b, a);
}

int main() {
  cin.tie(nullptr)->sync_with_stdio(false);
  int p;
  cin >> p;
  a[0] = 1;  // 记录 a 数组的位数
  a[1] = 2;  // 对 2 进行平方
  b[0] = 1;  // 记录 b 数组的位数
  b[1] = 1;  // 答案数组
  binpow(p);
  cout << (int)(log10(2) * p) + 1 << '\n';
  b[1] -= 1;  // 最后一位减 1
  for (int i = M; i >= 1; i--) {
    cout << b[i];
    if ((i - 1) % 50 == 0) {
      cout << '\n';
    }
  }
}

Tiền xử lý lũy thừa nhanh khi cố định cơ số

Khi cơ số \(a\) cố định, có thể dùng phân khối để tiền xử lý và trả lời mỗi truy vấn trong \(O(1)\). Thuật toán này còn gọi là “lũy thừa ánh sáng”. Quy trình:

  1. Chọn \(s\), tiền xử lý \(a^0,a^1,\cdots,a^{s-1}\) và \(a^0,a^s,\cdots,a^{\lfloor p/s\rfloor s}\) vào hai mảng;
  2. Với truy vấn \(a^b\), tách \(b=\lfloor b/s\rfloor s+(b\bmod s)\), khi đó \(a^b=a^{\lfloor b/s\rfloor s}\cdot a^{b\bmod s}\), trả lời \(O(1)\).

Nếu \(b \in [0,n]\), thường chọn \(s\approx \sqrt{n}\) hoặc gần lũy thừa của \(2\). Chọn \(\sqrt{n}\) tối ưu tiền xử lý \(O(\sqrt{n})\), còn chọn lũy thừa của \(2\) giúp đơn giản hóa bằng bit.

Đặc biệt với lũy thừa modulo, cơ số \(a\) cố định ngầm nghĩa mô-đun \(m\) cũng cố định. Do định lý Euler mở rộng, biên trên mũ là \(n = 2\varphi(m)\); với mô-đun nguyên tố \(p\), biên là \(n = p - 1\). Cả hai trường hợp đều tiền xử lý \(O(\sqrt{m})\).

Mã tham khảo 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
int a, mod, pow1[65536], pow2[65536];

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

Bài tập

Một phần nội dung trang này được dịch từ bài Бинарное возведение в степень và bản dịch tiếng Anh Binary Exponentiation. Bản tiếng Nga có giấy phép Public Domain + Leave a Link; bản tiếng Anh có giấy phép CC-BY-SA 4.0.

Last updated on this page:, Update history
Found an error? Want to help improve? Edit this page on GitHub!
Contributors to this page:OI-wiki
All content on this page is provided under the terms of the CC BY-SA 4.0 and SATA license, additional terms may apply