Berlekamp-Massey演算法學習筆記

Berlekamp-Massey演算法

很久之前就聽說過這個演算法，當時六校聯考的時候Day1T1是一道很有意思的遞推，神仙zzx不會做於是就拿BM演算法艹出了遞推式Orzzzzzzzzzzx

推薦一篇講的詳細的不能再詳細的博客

我就不詳細說了，只記一下自己感覺比較難理解的地方

設\(r(m)\)表示序列的遞推式且長度為\(m\)

\(f(r, i)\)表示\(\sum_{j = 1}^m r_j * a[i - j]\)

\(\delta(r, i)\)表示\(a[i] - f(r, i)\)

\(fail_i\)表示第\(i\)個遞推式出錯的位置

對於某一個位置\(i\)，如果我們求出的\(\delta(r, i) \not = 0\)，這時候我們需要構造一個遞推式\(r'(m')\)，滿足\(\forall j \in [m' + 1, i - 1] f(r', j) = 0\)且\(f(r, i) = \delta(r, i)\)

這樣我們令\(r = r + r'\)就得到新位置的遞推式了

\(r'\)可以這麼構造

設\(mul = \frac{\delta(r, i)}{\delta(r, fail_{cnt - 1})}\)

那麼\(r' = \{0, 0, 0 \dots, 0, mul, -mul * R_{cnt - 1} \}\)

\(0\)的個數為\(i - fail_{cnt - 1} - 1\)

至於為什麼這麼構造是對的，我思考了挺長時間，簡單的證明一下

首先對於\(\forall j \in [m' + 1, i - 1]\), \(\delta(r', j) = 0\)

仔細想了想，，發現自己並不會證。。如果哪位大佬會的話可以教教本蒟蒻

感性理解就是因為\(r\)在\([1, M]\)處滿足任意位置為\(0\)，然後右移一下還滿足？。。

至於為什麼\(f(r', i) = \delta(r, i)\)

可以這麼考慮，前\(i - fail_{cnt - 1} - 1\)個位置產生的貢獻為\(0\)

\(mul\)產生的貢獻為\(mul * a_{fail_{cnt - 1}}\)

\(-mul * R_{cnt - 1}\)產生的貢獻為\(-mul * (a[fail_{cnt - 1}] - \delta(r, fail_{cnt - 1]})\)

合併同類項後可以得到\(mul * \delta(r, fail_{cnt - 1}) = \delta(r, i)\)

代碼如下

#include<bits/stdc++.h>
using namespace std;
const int MAXN = 2005;
const double eps = 1e-8;
int cnt, fail[MAXN];
double val[MAXN], delta[MAXN];
vector <double> ans[MAXN];
int main() {
    int N; scanf("%d", &N);
    for (int i = 1; i <= N; i++) scanf("%lf", &val[i]);
    for (int i = 1; i <= N; i++) {
        double tmp = val[i];
        for (int j = 0; j < ans[cnt].size(); j++)
            tmp -= ans[cnt][j] * val[i - j - 1];
        delta[i] = tmp;
        if (fabs(tmp) <= eps) continue;
        fail[cnt] = i;
        if (cnt == 0) {
            ans[++cnt].resize(i);
            continue;
        }
        double mul = delta[i] / delta[fail[cnt - 1]];
        cnt++; ans[cnt].resize(i - fail[cnt - 2] - 1);
        ans[cnt].push_back(mul);
        for (int j = 0; j < ans[cnt - 2].size(); j++)
            ans[cnt].push_back(ans[cnt - 2][j] * -mul);
        if (ans[cnt].size() < ans[cnt - 1].size()) ans[cnt].resize(ans[cnt - 1].size());
        for (int j = 0; j < ans[cnt - 1].size(); j++)
            ans[cnt][j] += ans[cnt - 1][j];
    }
    for (int i = 0; i < ans[cnt].size(); i++)
        cout << ans[cnt][i] << ' ';
    return 0;
}