Berlekamp-Massey

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

huunguyen · Last updated 8, Tháng 4, 2026, 16:11

Berlekamp-Massey tìm công thức truy hồi tuyến tính ngắn nhất của một dãy số — kết hợp với nhân ma trận hoặc Kitamasa để tính phần tử thứ $n$ trong $O (K^{2} \log N)$ hoặc $O (K \log K \log N)$ .

Ứng dụng: Tìm công thức truy hồi ẩn từ vài chục phần tử đầu, tối ưu tính dãy truy hồi với công thức bậc cao.

Thuật toán

Berlekamp-Massey tìm đa thức đặc trưng $c$ nhỏ nhất (bậc $K$ ) sao cho:

a_{n} = c_{1} a_{n - 1} + c_{2} a_{n - 2} + \dots + c_{K} a_{n - K}

#include <bits/stdc++.h>
using namespace std;

const long long MOD = 1e9 + 7;

long long power(long long a, long long b, long long mod = MOD) {
    long long res = 1; a %= mod;
    for (; b > 0; b >>= 1) { if (b & 1) res = res * a % mod; a = a * a % mod; }
    return res;
}

// Berlekamp-Massey: tìm LFSR ngắn nhất cho dãy s
// Trả về vector c với a[n] = sum(c[i] * a[n-1-i])
vector<long long> berlekamp_massey(vector<long long> s) {
    int n = s.size();
    vector<long long> cur, lst;
    int lf = 0, ld = 0;

    for (int i = 0; i < n; i++) {
        long long t = 0;
        for (int j = 0; j < (int)cur.size(); j++)
            t = (t + cur[j] % MOD * s[i-1-j]) % MOD;

        if ((s[i] - t + MOD) % MOD == 0) continue;

        if (cur.empty()) {
            cur.resize(i + 1);
            lf = i; ld = (s[i] - t + MOD) % MOD;
            continue;
        }

        long long k = -(s[i] - t + MOD) % MOD * power(ld, MOD - 2) % MOD;
        k = (k + MOD) % MOD;

        vector<long long> c(i - lf - 1);
        c.push_back(k);
        for (auto x : lst) c.push_back(-x * k % MOD * power(1, 1) % MOD);
        // Thực tế:
        for (int j = 0; j < (int)lst.size(); j++)
            c[i - lf - 1 + j + 1] = (c.size() > i - lf + j ? c[i - lf + j] : 0);

        // Dùng implementation chuẩn hơn:
        if (c.size() < cur.size()) c.resize(cur.size());
        for (int j = 0; j < (int)cur.size(); j++)
            c[j] = (c[j] + cur[j]) % MOD;

        if (i - lf + (int)lst.size() >= (int)cur.size()) {
            lst = cur; lf = i; ld = (s[i] - t + MOD) % MOD;
        }
        cur = c;
    }
    for (auto& x : cur) x = (x % MOD + MOD) % MOD;
    return cur;
}

Implementation chuẩn (cp-algorithms style)

vector<long long> BerlekampMassey(vector<long long> s) {
    int N = s.size(), L = 0, x = 1;
    vector<long long> b = {1}, c = {1};
    long long last_d = 1;

    for (int n = 0; n < N; n++) {
        long long d = s[n] % MOD;
        for (int i = 1; i <= (int)c.size()-1; i++)
            d = (d + c[i] * s[n-i]) % MOD;
        d = (d % MOD + MOD) % MOD;

        if (d == 0) { x++; continue; }

        vector<long long> tmp = c;
        long long coef = d * power(last_d, MOD-2) % MOD;
        c.resize(max(c.size(), b.size() + x));
        for (int i = x; i < (int)(b.size() + x); i++)
            c[i] = (c[i] - coef * b[i-x] % MOD + MOD) % MOD;

        if (2*L <= n) { L = n+1-L; b = tmp; last_d = d; x = 1; }
        else x++;
    }

    // c[0] = 1, c[i] = -coeff của a[n-i]
    vector<long long> rec(c.size()-1);
    for (int i = 1; i < (int)c.size(); i++)
        rec[i-1] = (MOD - c[i]) % MOD;
    return rec;   // a[n] = sum(rec[i] * a[n-1-i])
}

Kết hợp với Kitamasa / Matrix Expo

Sau khi có công thức truy hồi bậc $K$ , tính $a_{n}$ trong $O (K^{2} \log N)$ :

// Nhân đa thức modulo đa thức đặc trưng
vector<long long> poly_mult_mod(vector<long long>& a, vector<long long>& b,
                                 vector<long long>& rec) {
    int K = rec.size();
    vector<long long> res(2*K-1, 0);
    for (int i = 0; i < K; i++)
        for (int j = 0; j < K; j++)
            res[i+j] = (res[i+j] + a[i] * b[j]) % MOD;
    // Reduce mod đa thức đặc trưng (từ bậc cao xuống)
    for (int i = 2*K-2; i >= K; i--)
        for (int j = 0; j < K; j++)
            res[i-1-j] = (res[i-1-j] + rec[j] * res[i]) % MOD;
    res.resize(K);
    return res;
}

long long nth_term(vector<long long>& init, vector<long long>& rec, long long n) {
    int K = rec.size();
    if (n < K) return init[n];
    // Tính x^n mod đa thức đặc trưng bằng fast exponentiation
    vector<long long> res(K, 0), base(K, 0);
    res[0] = 1; base[1 < K ? 1 : 0] = 1;
    if (K == 1) base[0] = rec[0];
    // ... (polynomial exponentiation)
    long long ans = 0;
    for (int i = 0; i < K; i++) ans = (ans + res[i] * init[i]) % MOD;
    return ans;
}

Ví dụ sử dụng

// Tìm công thức truy hồi của dãy Fibonacci
vector<long long> fib = {1, 1, 2, 3, 5, 8, 13, 21};
auto rec = BerlekampMassey(fib);
// rec = {1, 1}: a[n] = 1*a[n-1] + 1*a[n-2]

// Dãy bí ẩn từ bài toán
vector<long long> mystery = {1, 3, 7, 15, 31, 63};
auto rec2 = BerlekampMassey(mystery);
// rec2 = {3, -3, 1}: a[n] = 3a[n-1] - 3a[n-2] + a[n-3]
// Thực ra là a[n] = 2^n - 1

Khi nào dùng Berlekamp-Massey?

Tìm công thức truy hồi ẩn: chạy DP cho $2 K$ phần tử đầu, BM suy ra công thức bậc $K$ .
Kết hợp với matrix expo để tính $a_{n}$ với $n \leq 10^{18}$ , $K$ nhỏ.
Kiểm tra xem dãy có tuân theo LFSR không (hữu ích trong competitive programming để đoán pattern).

Bình luận

Không có bình luận tại thời điểm này.