Berlekamp Massey
(fps/berlekamp_massey.hpp)

• View this file on GitHub
• Last update: 2024-06-25 15:37:23+09:00
• Include: #include "fps/berlekamp_massey.hpp"

説明

$N$ 次線形回帰数列 $a$ の母関数を $A(x)$ とする。 $A(x)$ の先頭 $2N$ 項以上を入力として与えると、以下を満たす次数 $N$ のmonicな多項式 $C(x)$ を返す。時間計算量は $O(N^2)$ である。

\[[x^i] C(x) \times A(x) = 0 ( N \leq i )\]

Depends on

• modint/base.hpp

Required by

• Black Box Linear Algebra (matrix/black_box_linear_algebra.hpp)

Verified with

• test/math/Find_Linear_Recurrence.test.cpp
• test/matrix/Determinant_of_Sparse_Matrix.test.cpp
• test/yuki/yuki_1112.test.cpp

Code

#pragma once

#include <algorithm>
#include <vector>

#include "../modint/base.hpp"

namespace ebi {

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

}  // namespace ebi

#line 2 "fps/berlekamp_massey.hpp"

#include <algorithm>
#include <vector>

#line 2 "modint/base.hpp"

#include <concepts>
#include <iostream>
#include <utility>

namespace ebi {

template <class T>
concept Modint = requires(T a, T b) {
    a + b;
    a - b;
    a * b;
    a / b;
    a.inv();
    a.val();
    a.pow(std::declval<long long>());
    T::mod();
};

template <Modint mint> std::istream &operator>>(std::istream &os, mint &a) {
    long long x;
    os >> x;
    a = x;
    return os;
}

template <Modint mint>
std::ostream &operator<<(std::ostream &os, const mint &a) {
    return os << a.val();
}

}  // namespace ebi
#line 7 "fps/berlekamp_massey.hpp"

namespace ebi {

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

}  // namespace ebi

Back to top page