Berlekamp-Massey Algorithm
(math/berlekamp_massey.cpp)

• View this file on GitHub
• Last update: 2021-11-04 23:22:32+09:00

Description

Berlekamp-Massey のアルゴリズムは，与えられた数列を生成する最短の線形漸化式を求めるアルゴリズムである．

Operations

• vector<T> berlekamp_massey(vector<T> a)
  • $a$ を生成する最短の線形漸化式を1つ求める
  • 時間計算量: $O(n^2)$

Note

漸化式がわからない数列fの最初のいくつかを愚直に求めて，Berlekamp-Masseyで漸化式を得て，高速 Kitamasa 法で第k項を求めるという暴力的な方法がある．

auto coef = berlekamp_massey(f);
f.resize(coef.size());
reverse(coef.begin(), coef.end());
auto ans = kitamasa(f, coef, k);

Reference

• Berlekamp–Massey algorithm
• Linear Recurrence and Berlekamp-Massey Algorithm

Verified with

• test/yosupo/find_linear_recurrence.test.cpp

Code

#pragma once
#include <vector>

template <typename T>
std::vector<T> berlekamp_massey(const std::vector<T>& a) {
    int n = a.size();
    std::vector<T> B = {1}, C = {1};
    T b = 1;
    int L = 0, m = 1;
    for (int i = 0; i < n; ++i) {
        T d = a[i];
        for (int j = 1; j < (int) C.size(); ++j) {
            d += a[i-j] * C[j];
        }
        if (d == 0) {
            ++m;
        } else {
            auto tmp =  C;
            if (C.size() < m + B.size()) {
                C.resize(m + B.size());
            }
            T f = d / b;
            for (int j = 0; j < (int) B.size(); ++j) {
                C[m + j] -= f * B[j];
            }
            if (2 * L <= i) {
                L = i + 1 - L;
                b = d;
                B = tmp;
                m = 1;
            } else {
                ++m;
            }
        }
    }
    std::vector<T> ret(L);
    for (int i = 1; i <= L; ++i) {
        ret[i-1] = -C[i];
    }
    return ret;
}

#line 2 "math/berlekamp_massey.cpp"
#include <vector>

template <typename T>
std::vector<T> berlekamp_massey(const std::vector<T>& a) {
    int n = a.size();
    std::vector<T> B = {1}, C = {1};
    T b = 1;
    int L = 0, m = 1;
    for (int i = 0; i < n; ++i) {
        T d = a[i];
        for (int j = 1; j < (int) C.size(); ++j) {
            d += a[i-j] * C[j];
        }
        if (d == 0) {
            ++m;
        } else {
            auto tmp =  C;
            if (C.size() < m + B.size()) {
                C.resize(m + B.size());
            }
            T f = d / b;
            for (int j = 0; j < (int) B.size(); ++j) {
                C[m + j] -= f * B[j];
            }
            if (2 * L <= i) {
                L = i + 1 - L;
                b = d;
                B = tmp;
                m = 1;
            } else {
                ++m;
            }
        }
    }
    std::vector<T> ret(L);
    for (int i = 1; i <= L; ++i) {
        ret[i-1] = -C[i];
    }
    return ret;
}

Back to top page