浅谈线性递推与 BM 算法

前言

作为普通退役人 & 数学水平不足人，如果有错误欢迎来喷。

以下内容若无声明则均视为在模大质数意义下进行运算。

本文出现的题目题单

本文的一些示例代码

线性递推

这里我们将会讨论线性递推的一些性质于快速求其中一项的方法。

定义

对于一个数列 \{a\}，有 \forall n \ge k,a_n = \sum_{i = 1}^{k} a_{n - i} r_i 那么称这个数列满足 k 阶的常系数齐次线性递推关系，简称线性递推，\{a\} 为一个线性递推数列，有限数列 \{r\}为其的一个线性递推式。（以下称为定义 1 或者常见定义式）

更加形式化地，称对于数列 \{a\}，有 \forall n \ge k - 1,\sum_{i = 0}^{k - 1} a_{n - i} r_i = 0 则称 \{r\} 为 \{a\} 的一个线性递归式，若 r_0 = 0，称其为 线性递推式，称一个存在线性递归式的数列为线性递推数列，线性递归的阶数为 k - 1，可以发现这个定义就是上面的常见定义方式移项的结果。（以下称为定义 2）

下面如果没有说明，都会使用前者，即较为常见的定义。

线性递推的封闭性

对于线性递推数列 \{a\},\{b\}，有限数列 \{t_0,t_1,\cdots,t_{m - 1}\}，常数 c，有：

求线性递推数列其中一项

回到递推式：a_n = \sum_{i = 1}^{k} a_{n - i}r_i，如果把式子里每个 a_{n - i} 都以递推式展开，最后会得到以下形式：

\large{ a_n = \sum_{i = 0}^{k - 1} p_i a_i }

也就是说，如果把 a_0 \sim a_{k - 1} 当作一组基，我们可以用其线性表示数列每一项。

Lemma 1： 对于一组已知的 a_n = \sum_{i = 0}^{k - 1} p_i a_i，有 a_{n + m} = \sum_{i = 0}^{k - 1} p_i a_{i + m}。

考虑从已知的线性表示 a_m = \sum_{i = 0}^{k - 1} p_i a_i 和 a_n = \sum_{i = 0}^{k - 1} q_i a_i 得到 a_{n + m}。

首先将 a_{n + m} 用 Lemma 1 展开，得到：

\large \begin{aligned} a_{n + m} &= \sum_{i = 0}^{k - 1} p_i a_{i + n}\\ &= \sum_{i = 0}^{k - 1} p_i \sum_{j = 0}^{k - 1} q_j a_{i + j}\\ &= \sum_{i = 0}^{2k - 2} \sum_{j = 0}^{i} p_j q_{i - j} a_i \end{aligned}

于是我们得到了用 \mathcal{O} (k^2) 时间以 a_n,a_m 得到 a_{n + m} 的方法，只需要 a_{k} \sim a_{2k - 2} 的线性表示即可。

于是我们使用快速幂的思维，令一个线性表示为 a_1，另一个为 a_0，每次将前者从 a_n 变为 a_{2n} 然后以快速幂的方式就可以在 \mathcal{O} (k^2 \log n) 的时间内计算 a_n 了。

using i64 = long long;
#define add(x,y) ( (((x) += (y)) >= MOD) && ((x) -= MOD) )

int a[N],r[N];

inline void PolyMul(const int *p,const int *q,int *rs,int n) {
    static int t[N << 1];memset(t,0,(n + 1) << 3);
    for(int i = 0;i < n;++i) for(int j = 0;j < n;++j)
        add(t[i + j],(i64)p[i] * q[j] % MOD);
    for(int i = (n << 1) - 2;i >= n;--i) if(t[i])
        for(int j = 0;j < n;++j) if(r[j])
            add(t[i - j - 1],(i64)r[j] * t[i] % MOD);
    for(int i = 0;i < n;++i) rs[i] = t[i];
}

int LinearRecurrence(int n,int k) {
    static int tmp[N],res[N];
    tmp[1] = res[0] = 1;
    for(;n;n >>= 1) {
        if(n & 1) PolyMul(tmp,res,res,k);
        PolyMul(tmp,tmp,tmp,k);
    }
    int ans = 0;
    for(int i = 0;i < k;++i)
        add(ans,(i64)a[i] * res[i] % MOD);
    return ans;
}

然后我们可以发现，上式是一个卷积形式，使用多项式卷积和多项式取模可以优化到 \mathcal{O} (k \log k \log n)。

【模板】常系数齐次线性递推

从矩阵角度出发，使用 Cayley-Hamilton 定理会得到类似的结论，这里我们下面会稍微提到。

线性递推的应用

快速求出数列中的一项 【模板】常系数齐次线性递推。

加速线性递推形式的 DP 转移 [NOI2017] 泳池。

以及下面我们会提到的，搭配 Berlekamp-Massey 算法使用。

Berlekamp-Massey 算法

在这里我们将会讨论 Berlekamp-Massey 算法（以下简称 BM 算法）在求最短线性递推式以及其他一些方面的应用。

BM 算法流程

BM 算法会对于一个已知部分长 n 的数列 \{a\}，求出其每个前缀的最短线性递推式。

我们称 a_1 \sim a_i 得到的最短线性递推式为 r_i，特殊地，定义 r_0 为一个空数列。

考虑增量法，添加一个 a_i 时，我们已经得到的是 r_0 \sim r_{i - 1}，首先用 r_{i - 1} 计算一个可能的 a_i 的值记为 a'_i，令 \Delta = (a_i- a'_i) \mod p，如果 \Delta = 0，那么递推式到现在仍正确，令 r_i = r_{i - 1}。

否则考虑 r_{i - 1}，如果 r_{i - 1} 为空，r_{i} 即为\{0,0,0,\cdots,0\}，i 个 0。

非空时，考虑使用一个特殊构造的递推式来“修正” r_{i - 1} 以得到 r_i，那么如果有一个递推式 r_{fixed} 使得其计算出 a'_0 \sim a'_{i - 1} 均为 0，且计算出 a'_i = \Delta，就可以简单加减得到 r_i。在第一次的时候我们直接构造一个，后面我们发现可以用之前构造出来的进行一些移项即可再次利用，这时候我们为了线性递推式最短，只保留一个使得构造出来 r_{fixed} 最短的用于“修正”的递推式以保证移项后总得最短。

using i64 = long long;
#define add(x,y) ( (((x) += (y)) >= MOD) && ((x) -= MOD) )
#define red(x,y) ( (((x) -= (y)) < 0) && ((x) += MOD) )

inline i64 qpow(i64 a,i64 b = MOD - 2) {
    i64 res = 1;
    for(;b;b >>= 1) {
        if(b & 1) res = res * a % MOD;
        a = a * a % MOD;
    }
    return res;
}

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

EI 博客中介绍了一种以有理函数重建来求最短线性递推式的方法，可以到 EI 的 csdn 博客查看：最短线性递推式求解与有理函数重建-Entropy Increaser。

BM 求递推式所需的预处理长度

我们有了一个方法可以求出一个数列的最短线性递推式，现试证明对于一个未知线性递推数列需要最少多少项才可以求出一个保证准确的最短线性递推式。

令 l_i 表示 r_i 的阶数加一，或者上文提到的定义 2 中线性递推式的长度。

Lemma 2：如果 r_{i - 1} 不是 \{a_0 \sim a_i\} 的最短线性递推式，那么 l_i \ge \max (l_{i - 1},i + 1 - l_{i - 1})

考虑反证法，假设有 l_i \le i - l_{i - 1}。令 r_{i - 1} 为 \{p_0,p_1,\cdots,p_{l_{i - 1}}\}，r_{i} 为 \{q_0,q_1,\cdots,q_{l_{i}}\} 。

于是有：

\large\begin{aligned} - \sum_{j = 1}^{l_{i - 1}} p_j a_{i - j} &= \sum_{j = 1}^{l_{i - 1}} p_j \sum_{k = 1}^{l_i} q_k a_{i - j - k} (i - l_{i - 1} \ge l_i)\\ &= \sum_{k = 1}^{l_i} q_k \sum_{j = 1}^{l_{i - 1}} p_j a_{i - j - k}\\ &= -\sum_{k = 1}^{l_i} q_k a_{i - k} = a_i \end{aligned}

得到：r_{i - 1} 是 \{a_0 \sim a_i\} 的最短线性递推式，矛盾。

于是 Lemma 2 得证。

Lemma 3：对于一个 k 阶线性递推数列 \{a\}，求出 \{a_0 \sim a_{2k - 1}\} 的最短线性递推式就是原数列最短线性递推式。

如果有一个 t \ge 2k，满足 \{a_0 \sim a_{2k - 1}\} 的最短线性递推式不是 \{a_0 \sim a_t\} 的线性递推式，那么由 Lemma 2 线性递推式长度至少为 t - l_t + 1，显然大于 k，矛盾。

于是 Lemma 3 得证，我们只需要求出二倍于线性递推式长度即可。

BM 的应用

求线性递推式

在一些题目中，我们可以证明答案或我们要求的某个值满足线性递推关系，可以暴力预处理前面数项然后 BM 打出递推式。

[BJOI2019] 勘破神机

这道题常见的一个解法是第一类斯特林数推式子 + 扩域，但是我们可以用 BM 通过此题。

初步推导得到的式子都是 \sum_{i = 1}^{n} \binom{f_i}{k} 这样的形式，其中 \{f\} 为一个低阶的线性递推数列

考虑原式，前缀和相当于卷上 \{1,1,1,\cdots\}，可以摘掉。

剩下部分为：

\large\begin{aligned} \binom{f_i}{k} &= \frac{\prod_{j = 0}^{k - 1}f_i - j}{k!}\\ \end{aligned} 综上，$\sum_{i = 1}^{n} \binom{f_i}{k}$ 为一个线性递推数列，我们可以 $\mathcal{O} (k ^2)$ 暴力求出前面一定的项数然后丢进 BM 里。 BM 打表可知，最长的递推式长度不超过 $2507$，暴力求其二倍即数列前面约 $5000$ 项即可，时间复杂度 $\mathcal{O(k^2 + k \log k \log n)}$。 **例题** [CF717A Festival Organization](https://www.luogu.com.cn/problem/CF717A) [[THUSCH2017] 如果奇迹有颜色](https://www.luogu.com.cn/problem/P7454) #### 优化 DP 对于可以矩阵乘法优化的 DP，其一定可以用 BM 求解，下面证明，这里我们使用的是定义 2 的形式，方便推式子。 > Lemma 4：任意 $n$ 阶方阵 $M$ 的数列 $\{I,M,M^2,M^3,\cdots\}$ 都有不长于 $n$ 的线性递推式。 考虑一个 $n$ 阶方阵 $M$ 的数列 $\{I,M,M^2,M^3,\cdots\}$，我们首先求 $M$ 的特征多项式 $P$，由 Cayley-Hamilton 定理可知，特征多项式是 $M$ 的一个零化多项式，即 $P(M) = 0$ 且满足 $[x^n]P(x) = 1$，次数为 $n$ 次。 $$ \large\begin{aligned} P(M) &= 0\\ \sum_{i = 0}^{n} p_{n - i} M^i &= 0\\ \sum_{i = 0}^{n} p_{n - i} M^{i + j} &= 0\\ \sum_{i = 0}^{n} p_{i} M^{n + j - i} &= 0\\ \end{aligned} $$ 符合上面线性递推的定义式，于是其最短线性递推式不长于 $n$。 由线性递推的封闭性可知，原数列每个元素乘一个相同的 $n$ 行的列向量仍满足线性递推关系，于是可以矩阵快速幂优化的 DP 都可以套 BM。 **例题** [[SHOI2013] 超级跳马](https://www.luogu.com.cn/problem/P3990) [CF1117D Magic Gems ](https://www.luogu.com.cn/problem/CF1117D) [CF506E Mr. Kitayuta's Gift](https://www.luogu.com.cn/problem/CF506E) #### 稀疏线性方程组 有一个 $n \times n$ 的**满秩**矩阵 $A$ 和一个向量 $b$，求一个向量 $x$ 使得 $Ax = b$ ，$A$ 中元素个数为 $e$ ，剩下部分全空。 转化问题，相当于求 $x = A^{-1}b

根据 Lemma 4，我们可以求出 {b,Ab,A^2b,A^3b,\cdots} 的线性递推式记为 \{r_0,r_1,\cdots,r_{k - 1}\}，可知阶数不超 n。

回看定义 2 的式子可得：

\large \sum_{i = 0}^{k - 1} A^ibr_{k - i - 1} = 0\\ \sum_{i = 0}^{k - 1} A^{i - 1}br_{m - i - 1} = 0\\ A^{-1}b = -\frac{1}{r_{k - 1}} (\sum_{i = 0}^{k - 2} A^ibr_{m - i - 2})

现在只需求 {b,Ab,A^2b,A^3b,\cdots,A^{2n}b} 即可，考虑矩阵乘法时单独计算非 0 位置贡献，复杂度为 \mathcal{O} (ne)，总复杂度为 \mathcal{O} (n(n + e))。

帮助找规律

由于一些题使用 BM 复杂度过高，但是题目本身性质很强，可以考虑使用 BM 帮助找规律。

例题

冬至

在 BM+ 线性递推无法直接通过时，打表得到递推式满足 r_{k,i} = (k - i - 1)i!，于是把预处理的 k^2 复杂度摘掉了。

参考

毛啸 关于数列递归式的一些研究

钟子谦 两类递推数列的性质和应用

郭晓旭 杨宽 线性递推关系与矩阵乘法