orz myy and ZZQ

(其实在今天之前我从来没有听说过这玩意儿只是翻阅论文的时候有毛爷爷引进过一个没有中文资料的黑科技的印象)
其实也确实是一个黑科技……哪位神仙脑子一拍想出了这么神奇的算法…….(其实我更感兴趣怎么证明QAQ)

Berlekamp-Massey算法

也简记为BM算法(不是字符串匹配那个BM算法!)。用于在已知一组长度为$n$的数列项时求出此数列的最短线性递推式.显然已知的项数越多结果越准确。

线性递推式的定义

对一个数列$(x_1, x_2, …, x_{n – 1}, x_n)$,存在一组数列$(p_1, p_2, …, p_{m-1}, p_m)$,使得$\forall x_i(m + 1 \leq i \leq n)$,有$x_i = \sum_{j = 1}^m p_j * x_{i – j}$.则称其为数列$(x_1, x_2, … , x_{n – 1}, x_n)$的一种线性递推式.如果$m$是最小的,则称其为最短线性递推式。

Berlekamp-Massey算法运行过程

这没有办法,毕竟我只会实现

定义一个递推式代入第$i$项的值:

\begin{eqnarray}
g(i) =
\begin{cases}
0\ & 1 &\leq &i &\leq {m} \\
\sum_{j = 1}^m p_ix_{i – j} – x_i &m + 1 &\leq &i &\leq n
\end{cases}
\end{eqnarray}

现在开始求递推式.但是可以发现我们最开始是完全懵逼的,所以我们随意口胡一个递推式{}就好了。
由于不知道它是不是对的(如果是对的我建议你去买彩票),所以我们要一项一项地检验,于是一检验就出锅了。

形象地描绘一下现在的处境:我们有一个递推式$g$,它可以满足在前$i – 1$项值为$0$,但在第$i$项的值为$s$。
于是我们想要一个递推式$f$,他满足前$i – 1$项为$0$,第$i$项为$1$,那么我们要求的递推式就是$g – sf$。
如果我们曾经面临过这个问题,我们就会有一个$f$,它满足前$j – 1$项为$0$,第$j$项为$1$(不为$1$就除一除嘛)。发现没有?调整一下补几个$0$就行了.
如果有一堆$f$?显然是最近的最优。

具体实现

用具体的例子模拟一下:
一个数列$(1, 2, 4, 10, 24, 50, 116)$,求它的递推式。

首先设为$\{\}$,发现$1$,感觉我能怎么办我也很绝望,于是设置为$\{0\}$。
到了$2$,挂掉,曾今的经验告诉我们$x[1] = 1$,于是令$f = \{1\}$,递推式修改为$\{2\}$。
发现$4$,一切良好。
发现$10$,得到了$-2$(挂$*2$),曾今的经验告诉我们$x[2] = 2$,于是有$f = \{0, 0.5, 0\}$,递推式要加上$2f = \{0, 1, 0\}$,我们把递推式改为了$\{2, 1, 0\}$。
发现$24$,再次良好。
发现$50$,挂$*3$.得到$8$,曾今的经验又告诉我们:$-x[4] + 2 * x[3] = -2$,$f = \{0, 0.5, -1\}, 8f = \{0, 4, -8\}, g – 8f = \{2, -3, 8\}$。
发现$116$,挂$*4$,得到$-8$,刚刚挂的经历告诉我们:$-x[6] + 2 * x[5] + x[4] = 8$,有$f = \{-\frac{1}{8}, \frac{1}{4}, \frac{1}{8}, 0\}, g + 8f = \{1, -1, 7, 0\}$。
人做起来很累,但是电脑就很$nice$了啊。

这里有一份不保证正确的代码:

for (int i = 1; i <= n; ++i) {
    tmp = -a[i];
    for (int j = 0; j < ans[cnt].size(); ++j)
        tmp += ans[cnt][j] * a[i - j - 1];
    del[i] = tmp;
    if (fabs(tmp) < eps) continue;
    fail[cnt] = i;
    if (!cnt) {
        ans[++cnt].resize(i);
        continue;
    }
    mul = -tmp / del[fail[cnt - 1]];
    vector <double> &ls = ans[cnt - 1];
    vector <double> cur;
    cur.resize(i - fail[cnt - 1] - 1);
    cur.push_back(-mul);
    for (int j = 0; j < ls.size(); ++j) cur.push_back(mul * ls[j]);
    if (cur.size() < ans[cnt].size()) cur.resize(ans[cnt].size());
    for (int j = 0; j < ans[cnt].size(); ++j) cur[j] += ans[cnt][j];
    ans[++cnt] = cur;
}
分类: 文章

发表评论

电子邮件地址不会被公开。 必填项已用*标注

你是机器人吗? =。= *