玄學小記.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 \neq 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;
}