algorithm/berlekamp_massey.cpp

• View this file on GitHub
• Last update: 2021-12-28 21:38:32+09:00

Depends on

• utility/rep.cpp

Verified with

• test/berlekamp_massey.test.cpp

Code

#pragma once
#include <vector>
#include "../utility/rep.cpp"

template <class T> std::vector<T> berlekamp_massey(const std::vector<T>& a) {
    const int n = a.size();
    std::vector<T> c = {T(-1)}, c2 = {T(0)};
    T r2 = 1;
    int i2 = 0;
    for (const int i : rep(n)) {
        T r = 0;
        const int d = c.size();
        for (const int j : rep(d)) r += c[j] * a[i - j];
        if (r == T(0)) continue;
        T coeff = -r / r2;
        const int d2 = c2.size(), shift = i - i2 + 1;
        if (d2 + shift <= d) {
            for (const int j : rep(d2)) c[j + shift] += c2[j] * coeff;
        } else {
            std::vector<T> tmp = c;
            c.resize(d2 + shift);
            for (const int j : rep(d2)) c[j + shift] += c2[j] * coeff;
            c2 = std::move(tmp);
            i2 = i + 1;
            r2 = r;
        }
    }
    c.erase(c.begin());
    return c;
}

#line 2 "algorithm/berlekamp_massey.cpp"
#include <vector>
#line 2 "utility/rep.cpp"
#include <algorithm>

class Range {
    struct Iter {
        int itr;
        constexpr Iter(const int pos) noexcept : itr(pos) {}
        constexpr void operator++() noexcept { ++itr; }
        constexpr bool operator!=(const Iter& other) const noexcept { return itr != other.itr; }
        constexpr int operator*() const noexcept { return itr; }
    };
    const Iter first, last;

  public:
    explicit constexpr Range(const int first, const int last) noexcept : first(first), last(std::max(first, last)) {}
    constexpr Iter begin() const noexcept { return first; }
    constexpr Iter end() const noexcept { return last; }
};

constexpr Range rep(const int l, const int r) noexcept { return Range(l, r); }
constexpr Range rep(const int n) noexcept { return Range(0, n); }
#line 4 "algorithm/berlekamp_massey.cpp"

template <class T> std::vector<T> berlekamp_massey(const std::vector<T>& a) {
    const int n = a.size();
    std::vector<T> c = {T(-1)}, c2 = {T(0)};
    T r2 = 1;
    int i2 = 0;
    for (const int i : rep(n)) {
        T r = 0;
        const int d = c.size();
        for (const int j : rep(d)) r += c[j] * a[i - j];
        if (r == T(0)) continue;
        T coeff = -r / r2;
        const int d2 = c2.size(), shift = i - i2 + 1;
        if (d2 + shift <= d) {
            for (const int j : rep(d2)) c[j + shift] += c2[j] * coeff;
        } else {
            std::vector<T> tmp = c;
            c.resize(d2 + shift);
            for (const int j : rep(d2)) c[j + shift] += c2[j] * coeff;
            c2 = std::move(tmp);
            i2 = i + 1;
            r2 = r;
        }
    }
    c.erase(c.begin());
    return c;
}

Back to top page