闲话

？咋这么晚了

其实我写 BM 的原因是这样的
我不知道想啥就想到了线性递推了
然后想到线性递推就突然记起来 zky 代码里有个 BM 函数
当时看没咋注意后来才发现不对劲
然后又调出来细看了一遍
nnd，zky 这代码暴力算出了前 \(2k + 2\) 项之后直接扔进 BM 给整了个最简递推式出来
寄啊这 \(k\) 得升级到 1e5 级别啊

Berlekamp-Massey 算法

如何给你随便输入前几项的序列找到一个递推式？Berlekamp-Massey 算法提供了一种方案使得你可以在 \(O(n^2)\) 的复杂度内得到最简常齐次线性递推式，其中 \(n\) 是你输入的项数。

我们记我们有的数列是 \((a_1, a_2, \dots, a_n)\)，提取前 \(k\) 项得到的最简递推式是 \(R_k\)。令 \(R_0\) 为空。

假设现在我们已经得到了 \(R_1\sim R_{i - 1}\)，我们现在需要求出 \(R_i\)。

记 \(R_{i - 1}\) 带入 \(i\) 得到的值和真实的 \(a_i\) 的差值为 \(\Delta_{i - 1}\)，如果 \(\Delta_{i - 1} = 0\) 则显然有 \(R_i = R_i - 1\)。反之，如果 \(R_{i - 1}\) 为空则拓展 \(R_i\) 为 \(i\) 个 \(0\)，作为初始化。

我们接下来要考虑的就是 \(\Delta_{i - 1} \neq 0\) 且 \(R_{i - 1}\) 非空的情况，这时需要对 \(R_{i - 1}\) 调整得到 \(R_i\)。

我们发现，由于我们要的是线性递推式，我们只需要构造一个新的递推式使得这两个递推式按位加和的结果能拟合上 \(a_i\) 就行了。
能发现，我们需要构造的就是 \(R' = (r_1, r_2, \dots, r_k)\)，满足 \(\forall w\in (k, i), \ \sum_{j = 1}^k r_j\times a_{w - j} = 0\)，且带入 \(i\) 得到的值是 \(\Delta_{i - 1}\)。
这样 \(R_i = R_{i - 1} + R'\) 即可。

我们可以选择一个 \(w\in [0, i - 1)\)，假设 \(R_w = (r_1, r_2, \dots, r_t)\)，则可以构造

\[R' = (0, 0, \dots, 0, \frac{\Delta_{i - 1}}{\Delta_w}, -\frac{\Delta_{i - 1}}{\Delta_w}\times r_1, - \frac{\Delta_{i - 1}}{\Delta_w}\times r_2, \dots, - \frac{\Delta_{i - 1}}{\Delta_w}\times r_t)
\]

这里前面需要有 \(i - w - 2\) 个 \(0\)。

找一个最短的 \(R'\) 即可。由于我们只需要记录当前的 \(R_i\) 和最优的 \(R_w\)，总空间复杂度是 \(O(n)\) 的。时间复杂度显然是 \(O(n^2)\) 的。

板子直接套个常齐次线性递推即可。

code

inline int _R_Div(poly F, poly G, ll n) {
    using namespace polynomial::fast_number_theory_transform;
    int len = max(F.size(), G.size()); 
    int m = 1, o = 0;
    while (m < len) m <<= 1, ++ o;
    F.resize(1 << o + 1), G.resize(1 << o + 1);
    while (n > m) {
        ntt(F.data(), o + 1), ntt(G.data(), o + 1);
        for (register int i = 0; i < m * 2; ++ i) F[i] = 1ll * F[i] * G[i ^ 1] % mod;
        for (register int i = 0; i < m; ++ i) G[i] = 1ll * G[i << 1] * G[i << 1 | 1] % mod;
        intt(F.data(), o + 1);
        intt(G.data(), o);
        for (register int i = 0, j = n & 1; i < len; i ++, j += 2) F[i] = F[j];
        for (register int i = len, ed = 1 << o + 1; i < ed; ++ i) F[i] = G[i] = 0;
        n >>= 1;
    } G.resize(m); G = G.inv();
    int s = n; n = F.size() - 1, m = G.size() - 1;
    int a = max(0, s - m), b = min(s, n), ans = 0;
    for (register int i = a; i <= b; ++ i) ans = (ans + 1ll * F[i] * G[s - i]) % mod;
    return ans;
}
inline int linear_recur(ll n, int k, poly f, poly a) {
    poly F(k + 1); F[k] = 1; assert(a.size() >= k); a.resize(k);
    for (register int i = 1; i <= k; ++ i) assert(0 <= f[i] and f[i] < mod), F[k - i] = (f[i] == 0 ? 0 : mod - f[i]);
    F.rev(); f = (a * F).slice(a.degree());
    return _R_Div(f, F, n);
}
inline poly BM(int n, poly a) {
    a.f.emplace(a.f.begin(), 0);
    poly ans, lst; 
    i32 w = 0; ll delt = 0;
    for (register int i = 1; i <= n; ++ i) {
        ll tmp = 0;
        for (register int j = 0; j < ans.size(); ++ j) 
            tmp = (tmp + 1ll * a[i - j - 1] * ans[j]) % mod;
        if ((a[i] - tmp + mod) % mod == 0) continue;
        if (!w) {
            w = i, delt = a[i] - tmp + mod; if (delt >= mod) delt -= mod;
            for (register int j = i; j; -- j)
                ans.f.emplace_back(0);
            continue;
        } poly now = ans;
        ll mult = 1ll * (a[i] - tmp + mod) * qp(delt, mod - 2) % mod;
        if (ans.size() < lst.size() + i - w) 
            ans.resize(lst.size() + i - w);
        ans[i - w - 1] += mult; 
        if (ans[i - w - 1] >= mod) ans[i - w - 1] -= mod;
        for (register int j = 0; j < lst.size(); ++ j)
            ans[i - w + j] = (ans[i - w + j] - 1ll * mult * lst[j] % mod + mod) % mod;
        if (now.size() - i < lst.size() - w) 
            lst = now, w = i, delt = (a[i] - tmp + mod) % mod;
    } return ans;
} 
inline int recur_by_bm(ll n, int k, poly a) {
    poly f = BM(k, a);
    cout << f << '\n';
    if (f.size() == 1) return 1ll * a[0] * qp(f[0], n - 1) % mod;
    f.f.emplace(f.f.begin(), 0); 
    return linear_recur(n, f.degree(), f, a);
}