Berlekamp Massey (線形回帰数列の係数計算)
(library/polynomial/berlekamp_massey.hpp)

• View this file on GitHub
• Last update: 2024-01-30 20:59:32+09:00
• Include: #include "library/polynomial/berlekamp_massey.hpp"

Berlekamp Massey (線形回帰数列の係数の計算)

体 $F$ 上の数列 $s=(s_0,s_1,\ldots,s_{n-1})$ に対して、

\[s_i=\sum_{j=1}^L c_j\cdot s_{i-j}\quad (\forall i\geq L)\]

を満たす数列 $c=(c_1,\ldots,c_L)$ のうち長さ $L$ が最小のものを $\Theta(n ^ 2)$ 時間で計算する。

Reference

[1] J. Massey, “Shift-register synthesis and BCH decoding,” in IEEE Transactions on Information Theory, vol. 15, no. 1, pp. 122-127, January 1969, doi: 10.1109/TIT.1969.1054260.

Verified with

• test/src/polynomial/berlekamp_massey/find_linear_recurrence.test.cpp

Code

#ifndef SUISEN_BERLEKAMP_MASSEY
#define SUISEN_BERLEKAMP_MASSEY

#include <cassert>
#include <vector>

namespace suisen {
    /**
     * @brief Find linear recurrence in O(|s|^2) time
     * @tparam F Arbitrary field (operator +, -, *, /, +=, -=, *=, /= must be defined)
     * @param s Prefix of a linearly reccurent sequence
     * @return The vector of length L+1 c s.t. c_0=1 and s_i=Sum[j=1,L]c_i*s_{i-j} for all i>=L, where L is the minimum integer s.t. there exists such c of length L+1.
     */
    template <typename F>
    std::vector<F> find_linear_recuurence(const std::vector<F>& s) {
        std::vector<F> B{ 1 }, C{ 1 };
        B.reserve(s.size()), C.reserve(s.size());
        F b = 1;
        std::size_t L = 0;
        for (std::size_t N = 0, x = 1; N < s.size(); ++N) {
            F d = s[N];
            for (std::size_t i = 1; i <= L; ++i) d += C[i] * s[N - i];
            if (d == 0) {
                ++x;
            } else {
                F c = d / b;
                if (C.size() < B.size() + x) C.resize(B.size() + x);
                if (2 * L > N) {
                    for (std::size_t i = 0; i < B.size(); ++i) C[x + i] -= c * B[i];
                    ++x;
                } else {
                    std::vector<F> T = C;
                    for (std::size_t i = 0; i < B.size(); ++i) C[x + i] -= c * B[i];
                    L = N + 1 - L, B = std::move(T), b = d, x = 1;
                }
            }
        }
        C.resize(L + 1);
        for (std::size_t N = 1; N <= L; ++N) C[N] = -C[N];
        return C;
    }
} // namespace suisen


#endif // SUISEN_BERLEKAMP_MASSEY

#line 1 "library/polynomial/berlekamp_massey.hpp"



#include <cassert>
#include <vector>

namespace suisen {
    /**
     * @brief Find linear recurrence in O(|s|^2) time
     * @tparam F Arbitrary field (operator +, -, *, /, +=, -=, *=, /= must be defined)
     * @param s Prefix of a linearly reccurent sequence
     * @return The vector of length L+1 c s.t. c_0=1 and s_i=Sum[j=1,L]c_i*s_{i-j} for all i>=L, where L is the minimum integer s.t. there exists such c of length L+1.
     */
    template <typename F>
    std::vector<F> find_linear_recuurence(const std::vector<F>& s) {
        std::vector<F> B{ 1 }, C{ 1 };
        B.reserve(s.size()), C.reserve(s.size());
        F b = 1;
        std::size_t L = 0;
        for (std::size_t N = 0, x = 1; N < s.size(); ++N) {
            F d = s[N];
            for (std::size_t i = 1; i <= L; ++i) d += C[i] * s[N - i];
            if (d == 0) {
                ++x;
            } else {
                F c = d / b;
                if (C.size() < B.size() + x) C.resize(B.size() + x);
                if (2 * L > N) {
                    for (std::size_t i = 0; i < B.size(); ++i) C[x + i] -= c * B[i];
                    ++x;
                } else {
                    std::vector<F> T = C;
                    for (std::size_t i = 0; i < B.size(); ++i) C[x + i] -= c * B[i];
                    L = N + 1 - L, B = std::move(T), b = d, x = 1;
                }
            }
        }
        C.resize(L + 1);
        for (std::size_t N = 1; N <= L; ++N) C[N] = -C[N];
        return C;
    }
} // namespace suisen

Back to top page