[可能有用科技]线性递推与BM算法

前言

不会线性代数。

在某次模拟结束后看题解，“用BM算法求出递推式即可” 这句风轻云淡的话极大伤害了我这个数学弱菜。

但是起码当时我还是知道这里的 BM 说的一定不是 Boyer-Moore 字符串匹配，不过光凭BM算法这个关键字似乎只能搜到 Boyer-Moore，而加上递推之类的关键字才可以搜到 Berlekamp–Massey 算法。

线性递推

对于一个数列，其有 k 阶线性递推关系如下 ：

f_n = \sum_{i = 1}^{k} a_i f_{n - i}

或者形象地说，是由之前求出的 k 个数乘上一个系数再加起来求得的。

对于这个问题，最简单的做法就是根据式子，直接 \mathcal{O} (nk) 求就可以。

但是人不能满足于此，于是可以通过把系数填到矩阵里然后用矩阵乘法来加速这个过程。

转移矩阵 k 行 k 列，长这个样子 ：

\begin{bmatrix} & 1& & & & &\\ & & 1& & & &\\ & & & \ddots& & &\\ & & & & 1& &\\ a_k& a_{k - 1}& \cdots& a_3& a_2& a_1& \end{bmatrix}

以下记为矩阵 T 吧。

那这个矩阵有啥用呢？

首先我们有数列中一些项的值，将其写入列矩阵中 ：

V_t = \begin{bmatrix} f_t\\ f_{t + 1}\\ f_{t + 2}\\ \vdots\\ f_{t + k - 1} \end{bmatrix}

然后转移矩阵左乘这个矩阵 ：

\begin{bmatrix} & 1& & & & &\\ & & 1& & & &\\ & & & \ddots& & &\\ & & & & 1& &\\ a_k& a_{k - 1}& \cdots& a_3& a_2& a_1& \end{bmatrix} \times \begin{bmatrix} f_t\\ f_{t + 1}\\ f_{t + 2}\\ \vdots\\ f_{t + k - 1} \end{bmatrix} = \begin{bmatrix} f_{t + 1}\\ f_{t + 2}\\ f_{t + 3}\\ \vdots\\ f_{t + k} \end{bmatrix}

矩阵快速幂优化线性递推，得到了 \mathcal{O} (k^3 \log n) 的复杂度。

但是缺点也是很明显的，对于斐波那契数列这样阶数小的递推，这无疑是非常快的，但是如果 k 是一个较大的值，或者说只是 10^3 \sim 2 \times 10^3 就已经使得矩阵快速幂变得不够优秀。

最好的自然是能直接求出数列的通项，但是人毕竟不能指望一口吃成个胖子。

考场上让你硬做生成函数可能比生个孩子还难。(—— \mathrm{{\color{black}w}{\color{red}ind\_whisper}} 语）

所以我们仍然看这个矩阵乘法的做法。

矩阵乘法都求了什么信息？

V_{n - k + 1} = \begin{bmatrix} f_{n - k + 1}\\ f_{n - k + 2}\\ f_{n - k + 3}\\ \vdots\\ f_{n} \end{bmatrix}

前面的 f_{n - k + 1} \sim f_{n - 1} 的部分有用吗 ？没用。

求出来 T^n 那么一大块 k \times k 的矩阵有用吗 ？没用。

懂了，都是因为有没用的信息才变慢的，那么就考虑不要没有用的信息了。

多项式视角

首先直接返璞归真 ：

看我们最开始有什么 ： V_1。

然后是要求什么 ： f_n。

那么倒着推 ：

f_n = \sum_{i = 1}^{k} a_i f_{n - i}

那么 f_{n - 1},f_{n - 2},\dots,f_{n - k} 也可以如是用下标更小的数列中的数表示，最终回归到 ：

f_n = \sum_{i = 1}^{k} p_if_i

那么可以把数列里每个元素当成一个多项式。

将

\sum p_i f_i

表示为

\sum p_i x^i

然后再看原来的式子

f_n = \sum_{i = 1}^{k} a_i f_{n - i}

表示成多项式的话就是一个高次多项式分成了低次项，类似多项式取模。

尝试构造多项式 \lambda 使得 ：

x^t \equiv \sum_{i = 1}^{k} p_ix^{t - i}\pmod\lambda

然后推式子 ：

&x^t \equiv \sum_{i = 1}^{k} p_ix^{t - i}\pmod\lambda\\ &x^t - \sum_{i = 1}^{k} p_ix^{t - i} \equiv 0\pmod\lambda\\ &x^{t - k}(x^k - \sum_{i = 1}^{k} p_ix^{k - i}) \equiv 0 \pmod\lambda \end{aligned}

于是 ：

\lambda = x^k - \sum_{i = 1}^{k} p_ix^{k - i}

全程在模 \lambda 意义下进行即可。

复杂度 \mathcal{O}(k \log n \log k)

这个做法的好处是即使你不写 NTT 也可以达到 \mathcal{O} (k^2 \log n)，并且可以像矩阵乘法那样任意模数。

代码

多项式缺省源

常数大大大大大！！！！

inline int LinearRecurrence(const int *a,const int *f,int n,int k) {
    static int L[N],F[N],A[N],Q[N],R[N];
    repb(i,k - 1,0) L[i] = (MOD - f[k - i - 1]) % MOD;
    L[k] = 1,F[0] = 1,A[1] = 1;
    for(;n;n >>= 1) {
        if(n & 1) {
            PolyMul(F,A,F,k,k);
            PolyDiv(F,L,Q,R,k << 1,k + 1);
            repl(i,0,k) F[i] = R[i];
        }
        PolyMul(A,A,A,k,k);
        PolyDiv(A,L,Q,R,k << 1,k + 1);
        repl(i,0,k) A[i] = R[i];
    }
    ll ans = 0;
    repl(i,0,k) ans = (ans + (ll)F[i] * a[i] % MOD) % MOD;
    return ans;
}

int f[N],a[N];

int main() {
    init_IO();
    init_root();
    int n = read(),k = read();
    repl(i,0,k) f[i] = ((ll)read() + MOD) % MOD;
    repl(i,0,k) a[i] = ((ll)read() + MOD) % MOD;
    write(LinearRecurrence(a,f,n,k)),enter;
    end_IO();
    return 0;
}

于是就有些题非常毒瘤，需要

那我说的是个什么题呢？

[NOI2017] 泳池

模数 ： 998244353，一看就是要搞事情了。

首先恰好如何的概率不好算，那就把恰好为 k 的概率通过差分变成至少为 k 的概率减去至少为 k - 1 的概率。

然后想个办法把 DP 转移式搞出来 :

求递推式

现在问题变得比较奇怪，我们有一个数列并且 猜测(?) 证明(?) 无论如何其为一个满足线性递推关系的数列。

然后考虑把 最短的 能描述递推关系的递推式求出来。

那么我们可以乱搞啊！

首先瞎掰一个递推式，在加入一个数使得不满足目前的递推关系时吃后悔药把瞎掰的假递推式改吧改吧变成真的，重复这个过程。

现在定义目前的递推式有 k 项，为 ：

f_n = \sum_{i = 1}^{k} p_if_{n - i}

然后定义一次代入得到的值为 ：

\sum_{i = 1}^{k} p_if_{n - i} - f_n

可以知道如果递推式合法则每个代入得到的值均为 0。

假设现在就进行到 n - 1 位全部没事，就是在第 n 位代入的时候出现问题了 ：

\sum_{i = 1}^{k} p_if_{n - i} - f_n = s

出现了一个非 0 的数 s。

那么如果现在有一个新的递推式满足前 n - 1 项且第 n 项结果为 1 那么就可以把这个递推式乘上 s 然后减去就行了。

然后我们发现这个递推式求出来之后可以再利用，先移项再补 0 即可。

时间复杂度 \mathcal{O} (n^2)

代码

inline int LinearRecurrence(const int *a,const int *f,int n,int k) {
    static int L[N],F[N],A[N],Q[N],R[N];
    repb(i,k - 1,0) L[i] = (MOD - f[k - i - 1]) % MOD;
    L[k] = 1,F[0] = 1,A[1] = 1;
    for(;n;n >>= 1) {
        if(n & 1) {
            PolyMul(F,A,F,k,k);
            PolyDiv(F,L,Q,R,k << 1,k + 1);
            repl(i,0,k) F[i] = R[i];
        }
        PolyMul(A,A,A,k,k);
        PolyDiv(A,L,Q,R,k << 1,k + 1);
        repl(i,0,k) A[i] = R[i];
    }
    ll ans = 0;
    repl(i,0,k) ans = (ans + (ll)F[i] * a[i] % MOD) % MOD;
    return ans;
}

inline int BerlekampMassey(const int *f,int *p,int n) {
    static int lstf[N],curf[N],tmpf[N];
    int w = 0,k = 0,lstk = 0;ll d = 0;
    rep(i,1,n) {
        ll tmp = 0;
        repl(j,0,k)
            tmp = (tmp + (ll)f[i - 1 - j] * curf[j]) % MOD;
        if(!((f[i] - tmp) % MOD)) continue;
        if(!w) {
            w = i,d = (ll)f[i] - tmp;
            repb(j,i,1)
                curf[k++] = 0;
            continue;
        }
        repl(j,0,k) tmpf[j] = curf[j];
        int tmpk = k;
        ll mul = ((ll)f[i] - tmp) * qpow(d) % MOD;
        if(k < lstk + i - w) k = lstk + i - w;
        curf[i - w - 1] = ((ll)curf[i - w - 1] + mul) % MOD;
        repl(j,0,lstk)
            curf[i - w + j] = ((ll)curf[i - w + j] - mul * lstf[j]) % MOD;
        if(tmpk - i < k - w) {
            repl(j,0,tmpk) lstf[j] = tmpf[j];
            lstk = tmpk,w = i,d = (ll)f[i] - tmp;
        }
    }
    repl(i,0,k) p[i] = (curf[i] % MOD + MOD) % MOD;
    return k;
}

int f[N],a[N];

int main() {
    init_IO();
    init_root();
    int n = read(),m = read();
    rep(i,1,n) f[i] = read();
    int k = BerlekampMassey(f,a,n);
    repl(i,0,k) write(a[i]),space;
    enter;
    write(LinearRecurrence(f + 1,a,m,k)),enter;
    end_IO();
    return 0;
}

然后就有了一个问题 ： 为什么需要这个自动求递推式的玩意呢？

假如我有一个 n^2 的 DP 如下 ：

f_{i,j} = \sum_{k} f_{i - 1,k}\times v_{j,k}

然后这个 v 是个和 f 没关系的可以单独求的量。

不出意外的话这时候就出意外了，你要求一个 n 极大的 f_{n,1}，显然直接 DP 是不行的，那就 BM 求出 f_{n,1} 的递推式然后就可以 \mathcal{O}(k^2 \log n) 求 （不会有题卡 \mathcal{O} (m^2 + k^2 \log n) 吧） f_{n,1} 了。