Lũy thừa nhị phân (binary exponentiation): tính nhanh aⁿ mod m trong C++ và Python

Vì sao nhân lặp n lần bị quá thời gian khi n = 10¹⁸?

Giả sử đề bài yêu cầu tính phần dư của aⁿ khi chia cho m = 10⁹ + 7 với số mũ n = 10¹⁸.

Nếu bạn dùng vòng lặp for chạy từ 1 đến n, mỗi bước nhân thêm a vào kết quả:

• Độ phức tạp thời gian là O(n). Một máy tính thực hiện được khoảng 10⁸ phép tính mỗi giây, nên với n = 10¹⁸ vòng lặp cần 10¹⁰ giây (hơn 300 năm) mới chạy xong!
• Bạn cũng không thể dùng std::pow(a, n) của C++ hay a ** n của Python rồi mới chia lấy dư % m: hàm std::pow trả về số thực double (chỉ chính xác 15–17 chữ số), còn a ** n trong Python sẽ cố tạo một số nguyên dài hàng tỷ chữ số làm tràn RAM ngay lập tức.

Ta cần một thuật toán vừa chạy trong O(log n), vừa chia lấy dư % m ở từng phép nhân trung gian.

Hai cách hiểu nguyên lý của Lũy thừa nhị phân

Có hai góc nhìn tương đương dẫn tới thuật toán lũy thừa nhị phân:

Góc nhìn 1: Tách số mũ chẵn / lẻ (Chia để trị)

Thay vì giảm số mũ đi 1 (aⁿ = a × aⁿ⁻¹), ta chia đôi số mũ dựa trên tính chẵn lẻ của n:

• Nếu n == 0: a⁰ = 1.
• Nếu n chẵn: Tính một lần giá trị nửa = aⁿ/², rồi tự nhân với chính nó: aⁿ = (aⁿ/²)².
• Nếu n lẻ: Tính a⁽ⁿ⁻¹⁾/², bình phương lên rồi nhân thêm a: aⁿ = (a⁽ⁿ⁻¹⁾/²)² × a.

Nhờ tính aⁿ/² đúng một lần rồi bình phương, sau mỗi bước số mũ n giảm đi một nửa!

Góc nhìn 2: Tách theo biểu diễn nhị phân của số mũ n

Mọi số nguyên dương n đều phân tích được thành tổng các lũy thừa của 2 dựa vào các bit 1 trong hệ nhị phân (xem bài phép toán bit bitwise).

Ví dụ với n = 13, biểu diễn nhị phân là 1101₂ vì 13 = 8 + 4 + 0 + 1. Khi đó: 3¹³ = 3⁸⁺⁴⁺¹ = 3⁸ × 3⁴ × 3¹

Các lũy thừa 3¹ → 3² → 3⁴ → 3⁸... được sinh ra liên tiếp bằng cách tự bình phương giá trị ngay trước đó. Do đó, ta chỉ việc duyệt từng bit của n từ phải sang trái:

• Nếu bit cuối của n bằng 1 (n & 1), nhân lũy thừa hiện tại vào biến kết quả kq.
• Bình phương cơ số a = a * a cho bit tiếp theo, và dịch phải n >>= 1 (chia nguyên n cho 2).

Ví dụ từng bước: Tính 3¹³ mod 1000

Giá trị đúng của 3¹³ là 1.594.323, nên 3¹³ mod 1000 = 323. Vòng lặp duyệt bit tính ra 323 chỉ sau 4 bước (kq = 1, a = 3, n = 13 = 1101₂):

Không có thời điểm nào giá trị trung gian vượt quá 1000 × 1000 = 10⁶ vì sau mỗi phép nhân ta đều lấy % 1000 ngay!

Tính chất đồng dư của phép nhân và cảnh báo tràn số trong C++

Cơ sở toán học cho phép lấy phần dư % m ở mọi bước trung gian là tính chất đồng dư: (x × y) % m = ((x % m) × (y % m)) % m

Khi cài đặt bằng C++, có 2 quy tắc quan trọng để tránh tràn kiểu long long (≈ 9,22 × 10¹⁸):

1. Luôn gán a %= m ngay đầu hàm trước vòng lặp: Nếu a = 10¹² và m = 10⁹ + 7, quên dòng a %= m sẽ khiến phép nhân a * a = 10²⁴ ở bước đầu tiên bị tràn long long. Khi rút gọn a %= m trước, tích kq * a và a * a luôn nhỏ hơn m² ≈ 10¹⁸ < 2⁶³ - 1.
2. Khi m > 3 × 10⁹ (như m = 10¹⁵): Ngay cả khi a < m, tích a * a vẫn tràn long long. Trong C++ (GCC/Clang 64-bit), hãy ép kiểu sang số nguyên 128-bit trước khi nhân: (long long)((__int128)a * a % m).

Code Lũy thừa nhị phân trong C++17 và hàm pow(a, b, mod) của Python

Trong Python, ngoài việc tự viết vòng lặp để hiểu bản chất, hàm tích hợp sẵn pow(a, n, m) (truyền đủ 3 đối số) đã được viết bằng C bên dưới với thuật toán lũy thừa nhị phân tối ưu, chạy cực nhanh và tự động xử lý số nguyên lớn:

#include <iostream>

using namespace std;

// 1. Lũy thừa nhị phân bằng vòng lặp (theo các bit nhị phân của n)
long long binpow_mod(long long a, long long n, long long m) {
    a %= m; // Bắt buộc lấy dư trước để a * a không tràn long long khi a >= 3e9
    long long kq = 1 % m;
    while (n > 0) {
        if (n & 1) { // Nếu bit cuối của n bằng 1 (n lẻ)
            kq = (kq * a) % m;
        }
        a = (a * a) % m; // Bình phương cơ số cho bit tiếp theo
        n >>= 1;         // Dịch phải 1 bit (tương đương n /= 2)
    }
    return kq;
}

// 2. Lũy thừa nhị phân bằng đệ quy (tách số mũ chẵn / lẻ)
long long binpow_de_quy(long long a, long long n, long long m) {
    if (n == 0) return 1 % m;
    long long nua = binpow_de_quy(a, n / 2, m);
    long long kq = (nua * nua) % m;
    if (n % 2 == 1) {
        kq = (kq * (a % m)) % m;
    }
    return kq;
}

int main() {
    const long long MOD = 1000000007LL;

    // 1. Ví dụ trong bài: 3^13 mod 1000 = 323
    cout << binpow_mod(3, 13, 1000) << " " << binpow_de_quy(3, 13, 1000) << "\n";

    // 2. Số mũ cực lớn n = 10^18 chỉ mất ~60 bước lặp
    cout << binpow_mod(2, 1000000000000000000LL, MOD) << "\n";

    // 3. Cơ số lớn a = 10^12 (nhờ a %= m đầu hàm nên không bị tràn số khi nhân a * a)
    cout << binpow_mod(1000000000000LL, 1000000000000000000LL, MOD) << "\n";
    return 0;
}

def binpow_mod(a: int, n: int, m: int) -> int:
    """Lũy thừa nhị phân bằng vòng lặp trong O(log n)."""
    a %= m
    kq = 1 % m
    while n > 0:
        if n & 1:
            kq = (kq * a) % m
        a = (a * a) % m
        n >>= 1
    return kq


def binpow_de_quy(a: int, n: int, m: int) -> int:
    """Lũy thừa nhị phân bằng đệ quy."""
    if n == 0:
        return 1 % m
    nua = binpow_de_quy(a, n // 2, m)
    kq = (nua * nua) % m
    if n % 2 == 1:
        kq = (kq * (a % m)) % m
    return kq


if __name__ == "__main__":
    MOD = 1_000_000_007

    # 1. Ví dụ trong bài: 3^13 mod 1000 = 323
    print(binpow_mod(3, 13, 1000), binpow_de_quy(3, 13, 1000))

    # 2. Đối chiếu với hàm tích hợp sẵn pow(a, n, m) của Python
    n_lon = 10**18
    assert binpow_mod(2, n_lon, MOD) == pow(2, n_lon, MOD)
    print(pow(2, n_lon, MOD))

    a_lon = 10**12
    assert binpow_mod(a_lon, n_lon, MOD) == pow(a_lon, n_lon, MOD)
    print(pow(a_lon, n_lon, MOD))

Cả hai chương trình đều in ra:

• 323 323 (3¹³ mod 1000 bằng vòng lặp và đệ quy)
• 719476260 (2^(10¹⁸) mod 1000000007)
• 119177748 ((10¹²)^(10¹⁸) mod 1000000007).

Độ phức tạp thời gian và bộ nhớ

Vì số chữ số nhị phân của n là ⌊log₂ n⌋ + 1:

• Thời gian: O(log n) phép nhân. Với n = 10¹⁸, vì 2⁶⁰ > 10¹⁸, vòng lặp while (n > 0) chỉ chạy tối đa 60 bước!
• Bộ nhớ: O(1) với phiên bản dùng vòng lặp (và O(log n) bộ nhớ ngăn xếp nếu viết bằng đệ quy).

Ứng dụng: Tính nghịch đảo modulo và phép chia đồng dư (Định lý Fermat nhỏ)

Trong toán đồng dư, phép cộng, trừ và nhân đều tách lấy dư trực tiếp được, nhưng phép chia thì KHÔNG: (A / B) % m hoàn toàn khác với ((A % m) / (B % m)) % m!

Để thực hiện phép chia (A / B) mod p (ví dụ khi tính công thức tổ hợp C(n, k) = n! / (k! × (n - k)!) theo modulo p = 10⁹ + 7), ta phải nhân A với nghịch đảo modulo của B — tức số nguyên B⁻¹ sao cho (B × B⁻¹) % p == 1.

Theo Định lý Fermat nhỏ (Fermat's Little Theorem): Nếu p là một số nguyên tố và a không chia hết cho p (a % p != 0), thì: aᵖ⁻¹ ≡ 1 (mod p)

Tách aᵖ⁻¹ = a × aᵖ⁻² ≡ 1 (mod p), ta suy ra ngay công thức tính nghịch đảo modulo của a theo modulo nguyên tố p chính là một phép lũy thừa nhị phân: a⁻¹ ≡ aᵖ⁻² (mod p)

(Lưu ý: Nếu p không phải là số nguyên tố mà là hợp số, bạn không dùng được Định lý Fermat nhỏ mà phải dùng thuật toán Euclid mở rộng).

#include <iostream>

using namespace std;

long long binpow_mod(long long a, long long n, long long m) {
    a %= m;
    long long kq = 1 % m;
    while (n > 0) {
        if (n & 1) kq = (kq * a) % m;
        a = (a * a) % m;
        n >>= 1;
    }
    return kq;
}

// Nghịch đảo modulo của a theo số nguyên tố p (Định lý Fermat nhỏ: a^(p-2) mod p)
long long mod_inv_fermat(long long a, long long p) {
    return binpow_mod(a, p - 2, p);
}

// Tính tổ hợp C(n, k) % p bằng phép nhân với nghịch đảo modulo của mẫu số
long long to_hop_mod(int n, int k, long long p) {
    if (k < 0 || k > n) return 0;
    if (k > n - k) k = n - k;
    long long tu = 1, mau = 1;
    for (int i = 1; i <= k; ++i) {
        tu = (tu * (n - k + i)) % p;
        mau = (mau * i) % p;
    }
    return (tu * mod_inv_fermat(mau, p)) % p;
}

int main() {
    const long long MOD = 1000000007LL;

    // 1. Nghịch đảo của 5 theo modulo 10^9 + 7 là 400000003 (vì 5 * 400000003 % MOD == 1)
    long long inv5 = mod_inv_fermat(5, MOD);
    cout << inv5 << " " << ((5 * inv5) % MOD) << "\n";

    // 2. Tính C(10, 3) và C(100, 50) modulo 10^9 + 7
    cout << to_hop_mod(10, 3, MOD) << " " << to_hop_mod(100, 50, MOD) << "\n";
    return 0;
}

def mod_inv_fermat(a: int, p: int) -> int:
    """Nghịch đảo modulo theo Định lý Fermat nhỏ khi p là số nguyên tố."""
    return pow(a, p - 2, p)


def to_hop_mod(n: int, k: int, p: int) -> int:
    """Tính C(n, k) % p bằng cách nhân tử số với nghịch đảo modulo của mẫu số."""
    if k < 0 or k > n:
        return 0
    k = min(k, n - k)
    tu, mau = 1, 1
    for i in range(1, k + 1):
        tu = (tu * (n - k + i)) % p
        mau = (mau * i) % p
    return (tu * mod_inv_fermat(mau, p)) % p


if __name__ == "__main__":
    MOD = 1_000_000_007

    inv5 = mod_inv_fermat(5, MOD)
    print(inv5, (5 * inv5) % MOD)

    print(to_hop_mod(10, 3, MOD), to_hop_mod(100, 50, MOD))

Chương trình in ra:

• 400000003 1 (nghịch đảo của 5 theo modulo 10⁹ + 7 là 400000003, kiểm tra (5 × 400000003) % 1000000007 == 1).
• 120 538992043 (giá trị của C(10, 3) mod (10⁹ + 7) và C(100, 50) mod (10⁹ + 7)).

Lỗi hay gặp

Luyện ngay trên Big-O Coding

Thực hành ngay lũy thừa nhị phân và số học modulo qua các bài tập dưới đây trên Big-O Coding:

Câu hỏi thường gặp

Nguồn tham khảo

• cp-algorithms: Binary Exponentiation
• Python Documentation: Built-in function pow(base, exp, mod)
• CSES Competitive Programmer's Handbook — Chapter 21: Modular Arithmetic