玄学小记.6 ~ Berlekamp_Massey

？？？？？

怎么找一个数列的规律（线性递推）呢？当然就用BM啦！

估计这个东西我以后也遇不到几次……

为什么这个东西会出现在模拟赛里？？？

这个算法有什么用呢？
比如说有一道题，在 $m * n$ 的网格上搞一些事情，$m$ 非常小，$n$ 非常大。
显然是一个状压dp套矩阵快速幂~~裸题~~。
不过呢，虽然这个矩阵比较稀疏，但是边长有 $2000$ ……

这个时候就没法搞了 ……
不过可以猜测答案是一个关于 $n$ 的线性递推 ……
然后瞎搞搞发现这个线性递推只有 $1000$ 项，于是就可以做了 ……

BM就是一个用来求给定数列最短线性递推的算法。
~~但是为什么这个算法找到的是最短的啊？？？？~~

考虑一个长度为 $n$ 的数列 $S$，令 $Lambda$ 为数列 $S$ 的递推式，也就是说对于每个 $k geq L$ 都有
$$sum_{p = 0}^{L} Lambda_{p} S_{k - p} = 0$$

为了接下来好描述一些，考虑两个东西的生成函数 $Lambda(x) = sum_{x} Lambda_i x^i$ 与 $S(x) = sum_{x} S_i x^i$ 。上面的等式是个卷积，于是就有
$$Lambda(x)S(x) equiv M(x) pmod { x ^ n}$$
其中 $deg(M) < deg(Lambda)$ 。

如果我们找到了一个满足
$$C_{i - 1}(x)S(x) equiv D_{i - 1}(x) pmod {x^{i - 1}} 且 deg(D_{i - 1}) < deg(C_{i - 1})$$

的 $C_{i - 1}$ ，怎样求出一个满足

$$C_{i}(x)S(x) equiv D_{i}(x) pmod {x^{i}}且 deg(D_{i}) < deg(C_{i})$$
的 $C_{i}$ 呢？

显然

$$C_{i - 1}(x)S(x) equiv D_{i - 1}(x) + d x^{i - 1} pmod {x^i}$$

如果 $d = 0$ ，直接让 $C_i = C_{i - 1}$ ，否则我们需要想办法把这一项减掉。

如果我们找到了另一个 $C'_{i - 1}$，假设有
$$C'_{i - 1}(x)S(x) equiv D'_{i - 1}(x) + b x^{i - 1} pmod {x^i}$$
那么只需要令 $C_{i}(x) = C_{i - 1}(x) - frac{d}{b} C'_{i - 1}(x)$ ，这样 $dx^i$ 这一项就被消掉了。

怎么样得到一个 $C'_{i - 1}$ 呢？这里是BM算法的关键 —— 通过以前的信息去构造 $C'_{i - 1}$。
考虑满足
$$C_{j - 1}(x)S(x) equiv D_{j - 1}(x) + b x^{j - 1} pmod {mod x^j} 且 b eq 0$$
的 $C_{j - 1}$ ，也就是说它在递推时失配了。在等式两边同时乘上 $x^{i - j}$
$$x^{i - j}C_{j - 1}(x)S(x) equiv x^{i - j}D_{j - 1}(x) + b x^{i - 1}pmod {x^i}$$
于是我们就构造出了一个合法的 $C'_{i - 1}(x)$ !

因此
$$C_{i}(x) = C_{i - 1}(x) - frac{d}{b}x^{i-j}C_{j-1}(x)$$

注意到 $deg(C_i) - deg(C_{i-1}) leq 1$，因此 $j$ 只要取最近的一个就行了 ……

upd. 所以说要取的应该是$deg(C_{j}) - j$最小的$j$  …… 

递推的初始值是 $C_0 = 1, b = 1, i = 0, j = -1$ 。

好像非常好写的样子呢～

#include <bits/stdc++.h>
using namespace std;
const int L = 5000;
typedef vector <int> poly;
int p = 998244353;

void print(poly a)
{
    for (int i = 0; i < a.size(); ++ i) cerr << a[i] << " "; cerr << endl;
}
int powi(int a, int b)
{
    int c = 1;
    for (; b; b >>= 1, a = 1ll * a * a % p)
        if (b & 1) c = 1ll * c * a % p;
    return c;
}
poly operator + (poly a, poly b)
{
    poly c(max(a.size(), b.size()));
    for (int i = 0; i < a.size(); ++ i) c[i] = a[i];
    for (int i = 0; i < b.size(); ++ i) c[i] = (c[i] + b[i]) % p;
    return c;
}
poly operator - (poly a, poly b)
{
    poly c(max(a.size(), b.size()));
    for (int i = 0; i < a.size(); ++ i) c[i] = a[i];
    for (int i = 0; i < b.size(); ++ i) c[i] = (c[i] - b[i] + p) % p;
    return c;
}
poly operator << (poly a, int b)
{
    poly c(a.size() + b);
    for (int i = 0; i < a.size(); ++ i) c[i + b] = a[i];
    return c;
}
poly operator * (poly a, poly b)
{
    poly c(a.size() + b.size() - 1);
    for (int i = 0; i < a.size(); ++ i)
        for (int j = 0; j < b.size(); ++ j)
            c[i + j] = (c[i + j] + 1ll * a[i] * b[j]) % p;
    return c;
}
poly conv(int t)
{
    poly x(1); x[0] = t; return x;
}
poly Berlekamp_Massey(int S[], int len)
{
    poly Ci = conv(1), Cj = conv(1);
    int b = 1;
    for (int i = 0, j = -1; i < len; ++ i)
    {
        int d = 0;
        for (int j = 0; j < Ci.size(); ++ j)
            d = (1ll * Ci[j] * S[i - j] + d) % p;
        if (d)
        {
            poly tmp = Ci;
            Ci = Ci - ((conv(1ll * d * powi(b, p - 2) % p) * Cj) << (i - j));
            if ((int)Cj.size() - j > (int)tmp.size() - i)
                Cj = tmp, b = d, j = i;
        }
    }
    return Ci;
}
int X[L], len;
int main()
{
    cin >> len;
    for (int i = 0; i < len; ++ i) cin >> X[i];
    poly T = Berlekamp_Massey(X, len);
    for (int i = 0; i < T.size(); ++ i) cerr << T[i] << " "; cerr << endl;
}