Wiki › Thuật toán › Toán học › Berlekamp-Massey

Berlekamp-Massey

huunguyen Updated Tháng sáu 1, 2026

Berlekamp-Massey (BM) là thuật toán tìm công thức truy hồi tuyến tính (linear recurrence) ngắn nhất sinh ra một dãy số cho trước. Nói cách khác: cho vài chục phần tử đầu của một dãy $a_{0}, a_{1}, a_{2}, \dots$ , BM tìm bậc $K$ nhỏ nhất và các hệ số $r_{1}, \dots, r_{K}$ sao cho

$a_{n} = r_{1} a_{n - 1} + r_{2} a_{n - 2} + \dots + r_{K} a_{n - K} với mọi n \geq K .$

Khi đã có công thức truy hồi bậc $K$ , ta tính phần tử thứ $n$ (với $n$ rất lớn, vd $n \leq 10^{18}$ ) trong $O (K^{3} \log n)$ bằng nhân ma trận lũy thừa (matrix exponentiation), hoặc $O (K^{2} \log n)$ / $O (K \log K \log n)$ bằng Kitamasa. So với việc lặp $O (n)$ phần tử thì đây là bước nhảy vọt khi $n$ khổng lồ và $K$ nhỏ.

Điểm mạnh thực chiến: nhiều bài đếm/đại số cho ra một dãy mà ta không biết công thức, nhưng nếu sinh được $2 K$ phần tử đầu bằng vét cạn / DP, BM sẽ "đoán" ra công thức truy hồi, rồi matrix-expo nhảy tới phần tử thứ $n$ . BM chạy trong $O (N \cdot K)$ với $N$ là số phần tử đưa vào.

Ý tưởng / Trực giác

BM xử lý dãy tăng dần từng phần tử và duy trì một bộ hệ số ứng viên — gọi là đa thức đặc trưng $C (x) = 1 - c_{1} x - c_{2} x^{2} - \dots - c_{L} x^{L}$ bậc $L$ — sao cho công thức truy hồi tương ứng đúng cho toàn bộ tiền tố đã đọc.

Ý tưởng cốt lõi nằm ở khái niệm độ lệch (discrepancy) $d$ . Giả sử ta đang ở phần tử thứ $n$ và đang dùng truy hồi hiện tại có bậc $L$ . Ta tính giá trị mà truy hồi dự đoán cho $a_{n}$ , rồi lấy

$d = a_{n} - (c_{1} a_{n - 1} + c_{2} a_{n - 2} + \dots + c_{L} a_{n - L}) .$

Nếu $d = 0$ : truy hồi hiện tại vẫn đúng tới $a_{n}$ , không cần làm gì.
Nếu $d \neq 0$ : truy hồi sai ở vị trí này, phải sửa. Mấu chốt là ta đã lưu lại "lần sửa trước đó" — bộ hệ số $B (x)$ cùng độ lệch $d_{p r e v}$ tại bước mà nó lần đầu thất bại. Khi đó cộng thêm một bản dịch (shift) và co giãn của $B$ vào $C$ :

$C (x) - = \frac{d}{d_{p r e v}} \cdot x^{shift} \cdot B (x) .$

Tại sao phép cộng này lại "chữa" được đúng một độ lệch mà không phá các vị trí trước? Vì $B$ được chọn sao cho khi nhân với $x^{s h i f t}$ , nó tạo ra độ lệch $d_{p r e v}$ đúng tại vị trí $n$ và bằng 0 tại mọi vị trí trước đó. Nhân hệ số $d / d_{p r e v}$ làm độ lệch nó tạo ra đúng bằng $d$ , nên khi trừ đi, độ lệch tại $n$ bị triệt tiêu hoàn hảo.

Việc tăng bậc $L$ chỉ xảy ra khi $2 L \leq n$ . Đây là bất biến đảm bảo tính ngắn nhất: có một định lý (liên quan tới chặn dưới của độ phức tạp tuyến tính của LFSR) khẳng định rằng nếu một dãy tiền tố độ dài $n + 1$ không thể sinh bởi truy hồi bậc $\leq L$ thì bậc tối thiểu mới phải là $n + 1 - L$ . BM cập nhật $L \leftarrow n + 1 - L$ đúng theo công thức này, nên kết quả luôn là truy hồi bậc nhỏ nhất.

Trực giác tổng: BM giống như giải hệ phương trình tuyến tính tìm hệ số truy hồi, nhưng thực hiện gia tăng (incremental) — mỗi khi thêm một phương trình (một phần tử mới), nó chỉ chỉnh đúng phần cần chỉnh thay vì giải lại từ đầu, đạt $O (N K)$ thay vì $O (K^{3})$ của khử Gauss.

Ví dụ chạy tay

Lấy dãy $s = [1, 3, 7, 15, 31]$ (chính là $2^{n + 1} - 1$ ). Ta kỳ vọng BM tìm ra truy hồi bậc 2: $a_{n} = 3 a_{n - 1} - 2 a_{n - 2}$ .

Quy ước: $C$ là đa thức hiện tại (lưu dạng $[1, - c_{1}, - c_{2}, \dots]$ , tức $c_{0} = 1$ , $c_{i}$ là hệ số của $x^{i}$ ). $L$ là bậc truy hồi, $B$ là bản sao "lần sửa trước", $b$ là độ lệch lúc đó, $x$ là khoảng dịch.

Khởi tạo: $C = [1]$ , $L = 0$ , $B = [1]$ , $b = 1$ , $x = 1$ .

n=0: a=1
  d = s[0] = 1                  (chưa có hạng tử nào để trừ)
  d != 0, và 2L=0 <= n=0  -> tăng bậc
       C cũ = [1]; coef = d/b = 1/1 = 1
       C trở thành [1] - 1*x^1*[1] = [1, -1]
       L = 0+1-0 = 1 ;  B = [1] (C cũ) ; b = 1 ; x = 1
  Trạng thái:  C=[1,-1]  L=1
            chỉ số:   0    1
            ý nghĩa: a_n = 1*a_{n-1}

n=1: a=3
  dự đoán = c1*a0 = 1*1 = 1 ; d = 3 - 1 = 2
  d != 0, nhưng 2L=2 > n=1  -> KHÔNG tăng bậc, chỉ chỉnh C
       coef = d/b = 2/1 = 2
       C = [1,-1] - 2*x^1*[1] = [1, -1-2] = [1, -3]
       x++ -> x=2
  Trạng thái:  C=[1,-3]  L=1
            a_n = 3*a_{n-1}      (đúng tới n=1: 3 = 3*1 ✓)

n=2: a=7
  dự đoán = 3*a1 = 3*3 = 9 ; d = 7 - 9 = -2
  d != 0, và 2L=2 <= n=2  -> tăng bậc
       C cũ = [1,-3]; coef = d/b = -2/2 = -1   (b vẫn =1 từ n=0!)
       --- lưu ý: b=1 (độ lệch lần tăng bậc trước, ở n=0), x=2
       coef = d/b = -2/1 = -2
       C = [1,-3] - (-2)*x^2*[1] = [1, -3, +2] = [1,-3,2]
       L = 2+1-1 = 2 ; B = [1,-3] (C cũ) ; b = -2 ; x = 1
  Trạng thái:  C=[1,-3,2]  L=2
            a_n = 3*a_{n-1} - 2*a_{n-2}

  Kiểm tra: 3*7 - 2*3 = 21-6 = 15 = a3 ✓

n=3: a=15
  dự đoán = 3*a2 - 2*a1 = 21 - 6 = 15 ; d = 15-15 = 0  -> bỏ qua, x++

n=4: a=31
  dự đoán = 3*a3 - 2*a2 = 45 - 14 = 31 ; d = 0          -> bỏ qua

Kết quả: $C = [1, - 3, 2]$ , tức $a_{n} - 3 a_{n - 1} + 2 a_{n - 2} = 0$ , hay

$a_{n} = 3 a_{n - 1} - 2 a_{n - 2} .$

Chuyển sang dạng hệ số dương trả về: $r e c = [3, - 2]$ , đúng như kỳ vọng. (Trong cài đặt dưới đây mọi phép tính làm trên modulo, $- 2$ sẽ thành $MOD - 2$ .)

Cài đặt

Cài đặt BM trên trường $ℤ_{p}$ ( $p$ nguyên tố, dùng nghịch đảo modulo). Sau đó kết hợp với nhân ma trận để tính phần tử thứ $n$ .

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll MOD = 1e9 + 7;

ll power(ll a, ll b, ll mod = MOD) {
    ll res = 1; a %= mod; if (a < 0) a += mod;
    for (; b > 0; b >>= 1) {
        if (b & 1) res = res * a % mod;
        a = a * a % mod;
    }
    return res;
}
ll inv(ll a) { return power(a, MOD - 2); } // nghịch đảo: chỉ đúng khi MOD nguyên tố

// Trả về rec sao cho a[n] = sum_{i} rec[i] * a[n-1-i]
vector<ll> berlekampMassey(const vector<ll>& s) {
    vector<ll> C = {1}, B = {1};   // C: đa thức hiện tại; B: bản sao lần sửa trước
    int L = 0, m = 1;              // L: bậc truy hồi; m: khoảng dịch (shift)
    ll b = 1;                      // độ lệch tại lần tăng bậc gần nhất

    for (int n = 0; n < (int)s.size(); n++) {
        // tính độ lệch d = a[n] - (dự đoán theo truy hồi hiện tại)
        ll d = s[n] % MOD;
        for (int i = 1; i <= L; i++)
            d = (d + C[i] * s[n - i]) % MOD;
        d = (d % MOD + MOD) % MOD;

        if (d == 0) { m++; continue; }          // truy hồi vẫn đúng tới n

        if (2 * L <= n) {                        // cần tăng bậc
            vector<ll> T = C;                    // lưu C trước khi sửa
            ll coef = d % MOD * inv(b) % MOD;
            // C(x) -= coef * x^m * B(x)
            while ((int)C.size() < (int)B.size() + m) C.push_back(0);
            for (int i = 0; i < (int)B.size(); i++)
                C[i + m] = (C[i + m] - coef * B[i] % MOD + MOD) % MOD;
            L = n + 1 - L;                       // bậc tối thiểu mới
            B = T; b = d; m = 1;                 // cập nhật bản sao + độ lệch + reset dịch
        } else {                                 // chỉ chỉnh hệ số, không tăng bậc
            ll coef = d % MOD * inv(b) % MOD;
            while ((int)C.size() < (int)B.size() + m) C.push_back(0);
            for (int i = 0; i < (int)B.size(); i++)
                C[i + m] = (C[i + m] - coef * B[i] % MOD + MOD) % MOD;
            m++;
        }
    }

    // C = [1, -r1, -r2, ...]; trả rec[i] = -C[i+1]
    vector<ll> rec(C.size() - 1);
    for (int i = 1; i < (int)C.size(); i++)
        rec[i - 1] = (MOD - C[i]) % MOD;
    return rec;
}

// Nhân ma trận lũy thừa để tính a[n] từ truy hồi rec bậc K và K phần tử đầu
ll kth_term(const vector<ll>& rec, const vector<ll>& init, ll n) {
    int K = rec.size();
    if (n < (int)init.size()) return init[n] % MOD;
    // ma trận chuyển trạng thái K x K
    auto mul = [&](vector<vector<ll>>& A, vector<vector<ll>>& Bm) {
        vector<vector<ll>> R(K, vector<ll>(K, 0));
        for (int i = 0; i < K; i++)
            for (int k = 0; k < K; k++) if (A[i][k])
                for (int j = 0; j < K; j++)
                    R[i][j] = (R[i][j] + A[i][k] * Bm[k][j]) % MOD;
        return R;
    };
    vector<vector<ll>> M(K, vector<ll>(K, 0)), Res(K, vector<ll>(K, 0));
    for (int j = 0; j < K; j++) M[0][j] = rec[j];     // hàng đầu = hệ số truy hồi
    for (int i = 1; i < K; i++) M[i][i - 1] = 1;      // dịch trạng thái
    for (int i = 0; i < K; i++) Res[i][i] = 1;        // ma trận đơn vị
    ll e = n - (K - 1);
    for (; e > 0; e >>= 1) { if (e & 1) Res = mul(Res, M); M = mul(M, M); }
    // trạng thái ban đầu [a_{K-1}, ..., a_0]^T
    ll ans = 0;
    for (int j = 0; j < K; j++)
        ans = (ans + Res[0][j] * init[K - 1 - j]) % MOD;
    return ans;
}

Ví dụ dùng:

vector<ll> fib = {0, 1, 1, 2, 3, 5, 8, 13};
auto rec = berlekampMassey(fib);     // rec = {1, 1}: a_n = a_{n-1} + a_{n-2}
cout << kth_term(rec, {0, 1}, 100);  // F_100 mod 1e9+7

Độ phức tạp

Berlekamp-Massey: $O (N \cdot K)$ thời gian (chính xác hơn, tổng kích thước $C$ nhân số phần tử) với $N$ là số phần tử đưa vào và $K$ là bậc truy hồi tìm được, cộng $O (\log MOD)$ cho mỗi nghịch đảo modulo (có thể gộp lại còn vài lần). Mỗi phần tử $a_{n}$ phải duyệt qua $L \leq K$ hệ số để tính độ lệch, nên tổng là $O (N K)$ . Bộ nhớ $O (K)$ cho các đa thức $C, B$ cộng $O (N)$ lưu dãy đầu vào.
Nhân ma trận lũy thừa: $O (K^{3} \log n)$ thời gian — mỗi phép nhân hai ma trận $K \times K$ tốn $O (K^{3})$ , và lũy thừa cần $O (\log n)$ phép nhân. Bộ nhớ $O (K^{2})$ . Đây là lý do BM chỉ hiệu quả khi $K$ nhỏ (vài chục) còn $n$ rất lớn ( $10^{18}$ ).
Kitamasa (thay matrix-expo): $O (K^{2} \log n)$ , hoặc $O (K \log K \log n)$ nếu nhân đa thức bằng NTT — tốt hơn khi $K$ lớn.

Để BM tìm đúng truy hồi bậc $K$ , cần cấp ít nhất $2 K$ phần tử đầu (an toàn nên lấy nhiều hơn, vd $2 K + 5$ ); nếu không truy hồi có thể bị "khớp giả" với bậc thấp hơn thực tế.

⚠️ Lỗi thường gặp

Tràn số khi nhân: $a \cdot b$ với $a, b < 10^{9}$ vượt phạm vi int. Mọi tích phải ép long long và lấy modulo ngay: (a * b) % MOD. Quên % MOD trong vòng tính độ lệch là lỗi phổ biến nhất gây sai âm thầm.
Modulo âm: độ lệch $d$ và các hệ số sau phép trừ có thể âm. Luôn chuẩn hóa (x % MOD + MOD) % MOD. Đặc biệt khi trả về rec[i] = (MOD - C[i]) % MOD, nếu C[i] = 0 mà quên % MOD cuối thì ra MOD thay vì 0.
Nghịch đảo khi MOD không nguyên tố: inv(b) = power(b, MOD-2) chỉ đúng theo định lý Fermat khi MOD nguyên tố. Nếu bài cho modulo hợp số (vd $10^{9} + 9$ vẫn nguyên tố, nhưng $2^{32}$ thì không), power(b, MOD-2) cho kết quả sai — phải dùng nghịch đảo modulo qua Euclid mở rộng hoặc đổi cách tiếp cận. Ngoài ra nếu $b$ chia hết cho MOD thì không tồn tại nghịch đảo.
Cấp không đủ phần tử: đưa vào dưới $2 K$ phần tử khiến BM "đoán" ra truy hồi bậc thấp hơn thực, đúng với mẫu đã cho nhưng sai với phần tử kế tiếp. Luôn sinh dư phần tử và, nếu có thể, kiểm tra lại truy hồi tìm được trên vài phần tử cuối chưa dùng.
Sai chỉ số init khi matrix-expo: vector trạng thái ban đầu phải là $[a_{K - 1}, a_{K - 2}, \dots, a_{0}]^{T}$ theo đúng thứ tự ma trận; đảo thứ tự init hoặc nhầm số mũ $n - (K - 1)$ là bẫy off-by-one khiến lệch một phần tử.
Dãy không là truy hồi tuyến tính: BM luôn trả về một kết quả, kể cả khi dãy không hề tuân theo truy hồi tuyến tính cố định (vd dãy ngẫu nhiên) — khi đó bậc $K$ sẽ ≈ $N / 2$ và vô nghĩa. BM chỉ áp dụng khi bạn có cơ sở tin dãy là truy hồi tuyến tính (đếm cấu hình theo $n$ , lũy thừa ma trận chuyển trạng thái cố định...).

Biến thể / Mở rộng

Kitamasa thay nhân ma trận: tính $x^{n} mod C (x)$ trong vành đa thức để lấy $a_{n}$ , giảm từ $O (K^{3} \log n)$ xuống $O (K^{2} \log n)$ (hoặc $O (K \log K \log n)$ với NTT) — bắt buộc khi $K$ lớn.
BM trên GF(2): dùng cho dãy bit / thanh ghi dịch phản hồi tuyến tính (LFSR), cài bằng bitset, hữu ích trong mật mã và một số bài XOR.
Đoán công thức (guessing): mẹo thực chiến — viết brute-force sinh $30$ – $50$ phần tử đầu của đáp số theo $n$ , cho qua BM; nếu ra bậc nhỏ ổn định thì gần như chắc chắn đã tìm đúng truy hồi, rồi matrix-expo nhảy tới $n$ lớn. Liên hệ với nhân ma trận lũy thừa như công cụ tính phần tử thứ $n$ .

Bài tập luyện

Số Fibonacci (fibonum) — (Intermediate) Tính $F_{n}$ với $n \leq 10^{18}$ : bài "đích ngắm" kinh điển của truy hồi tuyến tính bậc 2, luyện đúng bước matrix-expo (hoặc BM tìm ra $a_{n} = a_{n - 1} + a_{n - 2}$ rồi nhảy tới $n$ ).
DZY và dãy Fibonacci (fibadd) — (Expert) Cập nhật đoạn cộng dãy Fibonacci và truy vấn tổng: khai thác cấu trúc truy hồi tuyến tính của Fibonacci để gộp cập nhật trên cây phân đoạn — luyện tư duy "Fibonacci như nghiệm của truy hồi tuyến tính".
Truy hồi tích lũy (prodrec) — (Expert) Truy hồi dạng tích $f_{x} = c^{2 x - 6} f_{x - 1} f_{x - 2} f_{x - 3}$ với $n \leq 10^{18}$ : lấy logarit (làm việc trên số mũ trong vành chỉ số) đưa về truy hồi tuyến tính bậc 3 cộng hạng tử đa thức, rồi matrix-expo — đúng tinh thần BM/Kitamasa khi $n$ khổng lồ.