[Asta 's Legacy]常系数齐次线性递推与 Berlekamp-Massey 算法

前言

本文是 “可能有用科技” 同名篇的优化版本。
写给冲刺 NOI2022 的机房的同志们。

因为笔者是线性代数渣渣，所以不会讲和线性代数相关的内容，即使讲了也大概率不会去证明。
想要看线性代数证明（特征多项式，Caylay-Hamilton 定理）的可以去看 yhx 神仙的线性代数分享。
笔者看了一点 zzq 神仙的集训队论文，但是由于太菜了没看懂，所以想深入学习的话笔者推荐 zzq 神仙的论文。
zzq 神仙的论文还有整式递推，过于高深.jpg
笔者看了一点毛爷爷的集训队论文，但是依然由于太菜了没看懂，所以想深入学习也可以去看看。
笔者只见过这个科技作为正解出在题里三次，分别是 [NOI2017] 泳池，[THUSCH2017] 如果奇迹有颜色，和一场 NOIP 前模拟赛的 T4。

常系数齐次线性递推

Linear recurrence with constant coefficients.
下文可能简称为线性递推。

“常系数”指系数均为常数。
“齐次”指没有常数项。

什么是线性递推

对于一个数列 \(\{f\}\)，满足 \(f_i = \sum_{j = 1}^{k} f_{i - j}a_j\) 的数列是一个常系数齐次线性递推的数列，此时我们称其为一个 \(k\) 阶的常系数齐次线性递推。

以下若无特殊声明，\(n\) 表示所求的某一项，\(k\) 表示线性递推阶数，数列下标从 \(0\) 开始，\(f_0 \sim f_{k - 1},a_1 \sim a_k\) 均为常数。

数列第 n 项的求法

法 1：按原式计算

既然已经给出了具体的计算方式，那么就可以 \(\mathcal{O}(nk)\) 地计算。

法 2:矩阵快速幂

通常遇到线性递推时都是 \(k\) 较小 \(n\) 较大的情况，那么为了优化，考虑让复杂度倾斜。
我们把已知的 \(k\) 项放进一个列矩阵 ：

\[\large F_0 = \begin{bmatrix} f_{k - 1}\\ f_{k - 2}\\ \vdots\\ f_0 \end{bmatrix} \]

我们希望得到的是：

\[\large F_{n - k + 1} = \begin{bmatrix} f_{n}\\ f_{n - 1}\\ \vdots\\ f_{n - k + 1} \end{bmatrix} \]

于是构造一个矩阵 \(T\) 使得 \(TF_{i - 1} = F_i\)：

\[\large T = \begin{bmatrix} &a_1\ &a_2 \ &a_3 &\cdots &a_k\\ &1 & & & & \\ & &1 & & & \\ & & &\ddots & & \\ & & & &1 &0 \end{bmatrix} \]

然后可知 \(T^n F_0 = F_n\)，根据矩阵乘法的结合律，先计算出 \(T^n\) 即可。
对于这部分，可以采取快速幂的思维，单次 \(k\) 阶方阵矩阵乘法复杂度为 \(k^3\)，总的时间复杂度为 \(k^3 \log n\)。

如果有多组询问该如何处理？
考虑预处理所有的 \(T^{2^i}\)，时间复杂度 \(k^3 \log n\)。
对于每次询问，只需要用一个列矩阵左乘即可，单次复杂度 \(k^2 \log n\)。

总体复杂度为 \(\mathcal{O} (k^3 \log n + Tk^2 \log n)\)。

代码没有，例题没有。

法 3：一个优化

观察原来的式子：

\[\large f_i = \sum_{j = 1}^{k} f_{i - j}a_j \]

可以发现，每个 \(f_n\) 都是由 \(f_{n - 1} \sim f_{n - k}\) 推出来的。
再用同样的方式把 \(f_{n - 1} \sim f_{n - k}\) 展开，最后式子必然形如：

\[\large \sum_{i = 0}^{k - 1} p_if_i \]

于是如果能快速地求出对于一个 \(f_n\)，将其表示出来的 \(p\) 的值，就可以快速计算了。

Lemma1：可加性

对于已知的

\[\large f_n = \sum_{i = 0}^{k - 1} p_i f_i \]

有

\[\large f_{n + x} = \sum_{i = 0}^{k - 1} p_i f_{i + x} \]

Lemma2: 由 \(f_n,f_m\) 到 \(f_{n + m}\)

现在已知：

\[\large{ f_n = \sum_{i = 0}^{k - 1} x_i f_i\\ f_m = \sum_{i = 0}^{k - 1} y_i f_i } \]

首先根据 Lemma 1：

\[\large f_{n + m} = \sum_{i = 0}^{k - 1} y_i f_{n + i} \]

然后可以把每一项 \(f_{n + i}\) 都展开：

\[\large{ f_{n + m} = \sum_{i = 0}^{k - 1} y_i \sum_{j = 0}^{k - 1} x_j f_{i + j} }\]

整理一下：

\[\large f_{n + m} = \sum_{i = 0}^{2k - 2} \sum_{j = 0} ^ {i} x_j y_{i-j} f_i \]

这个玩意是 \(\mathcal{O}(k^2)\) 的。
那么只需要已知 \(f_0 \sim f_{2k - 2}\) 就行了。
依然是按照快速幂的思想，求每个 \(f_{2^i}\) 的表示，最后就得到了 \(f_n\)。

时间复杂度 \(\mathcal{O} (k^2 \log n)\)，比起矩阵快速幂优化了一个 \(k\)。

代码：
有野蛮的宏定义，但是我没心思去改了。

其中 A[] 为系数。

#define add(x,y) ( (((x) += (y)) >= MOD) && ((x) -= MOD) )
int f[N],A[N];

inline void PolyMul(const int *a,const int *b,int *c,int n) {
	static int t[N];memset(t,0,(n + 1) << 3);
	repl(i,0,n) repl(j,0,n)
		add(t[i + j],(i64)a[i] * b[j] % MOD);
	for(int i = (n << 1);i >= n;t[i] = 0,--i) rep(j,1,n)
		add(t[i - j],(i64)A[j] * t[i] % MOD);
	repl(i,0,n) c[i] = t[i];
}

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

例题 [TJOI2010]天气预报。

目前最优解 50ms。

[NOI2017] 泳池

法4: 能不能再给力一点？

可以发现现在我们的复杂度瓶颈在于 \(\mathcal{O}(k^2)\) 的那一堆有点不知所云的 Lemma 2，把那个换成多项式乘法和多项式取模即可优化到 \(\mathcal{O} (k \log k \log n)\)

笔者写过但是因为码风过于丑陋就只挂个提交记录吧。

Link

例题没有。（怎么有人敢出只放 \(k \log k \log n\) 过的屑题啊）

二阶线性递推的特征方程法

对于一个二阶线性递推形如：

\[\large f_i = A f_{i - 1} + B f_{i - 2} \]

有特征方程：

\[x^2 - Ax - B = 0 \]

解出来两个根 \(x_0,x_1\)。
如果两根不同，有：

\[\large f_n = C x_0^n + D x_1^n \]

如果两根相同，有：

\[\large f_n = C x^n + Dn x^n \]

其中 \(C,D\) 为常数，代几个数进去算一下就行了。
这个方法最大的问题是会出现根式，那就需要根式里的常数在该题模数意义下具有二次剩余，这个可以考虑暴力找，最后多组询问靠光速幂解决即可。

例题 P5110 块速递推

为什么偏偏只说二阶的呢？因为更高阶做不到了，根据阿贝尔-鲁菲尼定理，五次及更高次的方程没有一般的根式解。

常系数齐次线性递推数列的前缀和

对于线性递推

\[\large f_i = \sum_{j = 1}^{k} f_{i - j}a_j \]

定义：

\[\large S_n = \sum_{i = 0}^{n} f_i \]

比较简单地，构造一个 \(k + 1\) 阶的转移矩阵就能解决这个问题了，复杂度 \(\mathcal{O} (k^3 \log n)\)。

（其实可以暴力算前面一部分前缀和然后丢进 BM 求递推式，可以证明是对的）

但是我不想写 \(k^2 \log\) 了，那就像 \(\mathrm{{\color{balck}{w}}{\color{red}{ind\_whisper}}}\) 一样挂个博客链接吧：https://www.luogu.com.cn/blog/lmpp/solution-p2461

例题就是这个 P2641，但是出题人害搁那惦记他那对合数取模呢，导致我的 BM 爆炸了。

Berlekamp-Massey 算法

很可惜的一点，通常我们是见不到把“我是常系数齐次线性递推哦”写在脸上的板子题的（要是敢出那出题人不被乱棍打死？），但是有一个十分强劲的算法，可以求一个数列的可能的最短递推式，所以可以在能证明其为线性递推的时候使用，当然不能证明也可以套一个上去骗分看看，就和不知道是不是多项式就套拉格朗日插值并且尝试多插几项一个道理。

这个算法是求模意义下的递推式的，所以如果遇到不当人的出题人用合数当模数那还是算了吧。

首先我们来骗来偷袭，最开始的递推式是空的，很快啊，这个递推式就要有东西了。
然后考虑如果到了第 \(i\) 个数的时候，现在我们有一个对于前 \(i - 1\) 个数都成立的递推式，用这个求出的第 \(i\) 项记为 \(f'_i\) ，然后实际的是 \(f_i\)，令 \(\Delta = (f'_i - f_i) \bmod p\)

如果 \(\Delta = 0\)，那现在的递推式就是对的。

否则继续考虑：
如果先前的递推式是空的，那现在就是 \(\{0,0,0,\cdots,0\}\) 了。
不然就要修正这个不优秀的递推式了。
那么如果现在有一个非常棒的递推式满足前 \(i - 1\) 项全 \(0\) 且第 \(i\) 项结果为 \(1\) 那么就可以把这个递推式乘上 \(\Delta\) 然后原式减去这个就行了。

那么我们考虑上一次 \(\Delta \ne 0\) 的时候，比如说是到了第 \(j\) 位的时候吧，那时留下了一个递推式恰好满足前 \(j - 1\) 为 \(0\) 然后 \(j\) 项为 \(\Delta_j\)，这时候除以一个 \(\Delta_j\) 然后稍微把这个式子移位一下就是我们要的一个递推式了，就叫它修正递推式吧。

然后我们只留着使得产生一个 “修正递推式” 最短的先前失配的递推式，这样就空间 \(\mathcal{O}(n)\)，时间 \(\mathcal{O} (n^2)\) 地解决这个问题了。

洛谷模板

代码：
码风依然非常的野蛮。

#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;
	rep(i,1,n) {
		i64 tmp = 0;
		repl(j,0,k)
			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;
			rep(j,1,i) curf[k++] = 0;
			continue;
		}
		repl(j,0,k) 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);
		repl(j,0,lstk)
			red(curf[i - w + j],(i64)lstf[j] * mul % MOD);
		if(tmpk - i < k - w) {
			repl(j,0,tmpk) lstf[j] = tmpf[j];
			lstk = tmpk,w = i,d = ((i64)f[i] - tmp + MOD) % MOD;
		}
	}
	repl(i,0,k) p[i] = ((i64)curf[i] + MOD) % MOD;
	return k;
}

int f[N],A[N];

inline void PolyMul(const int *a,const int *b,int *c,int n) {
	static int t[N];memset(t,0,(n + 1) << 3);
	repl(i,0,n) repl(j,0,n)
		add(t[i + j],(i64)a[i] * b[j] % MOD);
	for(int i = (n << 1);i >= n;t[i] = 0,--i) rep(j,1,n)
		add(t[i - j],(i64)A[j] * t[i] % MOD);
	repl(i,0,n) c[i] = t[i];
}

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

int main() {
	InitIO();
	int n = readI(),m = readI();
	rep(i,1,n) f[i] = readI();
	int k = BerlekampMassey(f,A + 1,n);
	rep(i,1,k) writeI(A[i]),space;
	enter;
	writeI(LinearRecurrence(m,k)),enter;
	enter;
	EndIO();
	return 0;
}

然后就有了一个问题 ： 为什么需要这个自动求递推式的玩意呢？
假如我有一个 \(n^2\) 的 DP 如下 ：

\[\large 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)\) 求了。

BM 还能做什么呢？
当你看到一个矩阵优化的题但是用矩阵快速幂复杂度高了那就用 BM 打出递推式吧。
因为线性递推的那个线性代数的做法就是理解为求 \(x^n\) 对转移矩阵特征多项式取模的结果，那我这么一个固定的转移矩阵也一定是可以这么求的。
例如 THUPC2022决赛 T4。

例题：

[SHOI2013]超级跳马
暴力 DP + BM，目前最优解。

[THUSCH2017] 如果奇迹有颜色
先 Polya 暴力 DP + BM 打表，然后有了表就线性递推即可。

[BJOI2019]勘破神机
然而正解不是 BM。
所以只放 \(k \log k \log n\) 的线性递推过

CF717A Festival Organization
上一题的子集，\(k^2 \log\) 可过。

CF506E Mr. Kitayuta's Gift
自动机上 DP 和矩阵 DP 狗都不写.jpg

[MtOI2019] 幻想乡数学竞赛
这是个数学题，但是 BM 可以打出线性递推式就是了。