常系数齐次线性递推优化
hly1204
·
2024-12-14 20:46:05
·
个人记录
看到有 UOJ 群友问常系数齐次线性递推怎么优化，时隔四年我来重置一下以前的文章（当然大部分还是拷贝的以前写的一些东西，好像没什么全新的东西）。
常系数齐次线性递推数列
我们一般称呼数列 \left(a_n\right)_{n\geq 0} 为 常系数齐次线性递推数列，当其满足下面的形式
a _ n=\sum _ {j=1}^d c _ j a _ {n-j}, \qquad (n \geq d)
其中 c_j 为常数并且不全为零，这时候我们称其为 d 阶的线性递推。给出 c _ 1, \dots , c _ d 和初值 a _ 0, \dots, a _ {d-1} ，我们的目标是在 k\gg d 时快速计算 a_k 。
此处我们定义 \Gamma(x):=x^d-\sum _ {j=0}^{d-1}c _ {d-j}x^j 为数列 \left(a_n\right)_{n\geq 0} 的 特征多项式 。系数我们考虑在 \mathbb{C} 上，当然本题中 \mathbb{F} _ {998244353} 上也是可以 FFT 的所以差不多。
友矩阵
对于首一多项式 \Gamma(x) 的友矩阵定义为
C_\Gamma:=
\begin{bmatrix}
&&&c _ d \\
1&&&c _ {d-1} \\
&\ddots &&\vdots \\
&&1&c _ 1
\end{bmatrix}
并且我们定义 b(x):=\sum _ {j=0}^{d-1}b _ jx^j 和
B _ b:=\begin{bmatrix}b _ 0 & b _ 1 & \cdots & b _ {d-1}\end{bmatrix}^{\intercal}
然后我们观察到将 C_{\Gamma} 应用在向量上的效果是
\underbrace{\begin{bmatrix}
&&&c _ d \\
1&&&c _ {d-1} \\
&\ddots &&\vdots \\
&&1&c _ 1
\end{bmatrix}} _ {C _ \Gamma}
\underbrace{\begin{bmatrix}
b _ 0 \\
b _ 1 \\
\vdots \\
b _ {d-1}
\end{bmatrix}} _ {B _ b}=
\underbrace{\begin{bmatrix}
c _ d b _ {d-1} \\
b _ 0 + c _ {d-1} b _ {d-1} \\
\vdots \\
b _ {d-2} + c _ 1 b _ {d-1}
\end{bmatrix}} _ {B _ {xb\bmod{\Gamma}}}
那么其自身的幂我们也可以定义了：
\begin{aligned}
C _ \Gamma &= \begin{bmatrix}B _ {x\bmod{\Gamma}} & B _ {x^2\bmod{\Gamma}} & \cdots & B _ {x^d\bmod{\Gamma}}\end{bmatrix}, \\
\left(C _ \Gamma\right)^2 &= \begin{bmatrix}B _ {x^2\bmod{\Gamma}} & B _ {x^3\bmod{\Gamma}} & \cdots & B _ {x^{d+1}\bmod{\Gamma}}\end{bmatrix}, \\
\vdots \\
\left(C _ \Gamma\right)^k&=\begin{bmatrix}B _ {x^k\bmod{\Gamma}} & B _ {x^{k+1}\bmod{\Gamma}} & \cdots & B _ {x^{k+d}\bmod{\Gamma}}\end{bmatrix}
\end{aligned}
降低阶
因为上面的线性递推是 d 阶的，d 有可能比较大，如果我们使用矩阵来描述，那么对于一个列向量而言，其就是 1 阶的。我们有
\begin{bmatrix}
a _ {k} \\
a _ {k+1} \\
\vdots \\
a _ {k+d-1}
\end{bmatrix}=\underbrace{\begin{bmatrix}
&1&& \\
&&\ddots & \\
&&&1 \\
c _ d&c _ {d-1}&\cdots &c _ 1
\end{bmatrix}^k} _ {\left(\left(C _ \Gamma\right)^{\intercal}\right)^k=\left(\left(C _ \Gamma\right)^{k}\right)^{\intercal}}
\begin{bmatrix}
a _ 0 \\
a _ {1} \\
\vdots \\
a _ {d-1}
\end{bmatrix}
那么使用矩阵快速幂就可以得到一个 O\left(\mathsf{MM}(d)\log k\right) 的算法，其中 \mathsf{MM}(d) 为计算两个 d\times d 的矩阵乘法的时间。
Fiduccia 算法
实际上我们在上面就可以观察到 \left(\left(C _ \Gamma\right)^{k}\right)^{\intercal} 的第一行也就是 \left(C _ \Gamma\right)^{k} 的第一列即 x^k\bmod{\Gamma(x)} 的系数，那么
a _ k=\left\langle x^k\bmod{\Gamma(x)},\begin{bmatrix} a _ 0 & \cdots & a _ {d-1}\end{bmatrix}\right\rangle
也就是两个向量的内积，这在上面矩阵描述中我们已经明确了。于是可以在 \mathbb{C}\left\lbrack x\right\rbrack /\left(\Gamma\right) 上计算快速幂。
形式 Laurent 级数
如果我们允许有限项的 x^{\lt 0} 的系数非零，那么我们就得到了形式 Laurent 级数
\mathbb{C}\left(\left(x\right)\right) := \left\lbrace \sum _ {j\geq N}a _ j x^j : N\in\mathbb{Z},a _ j \in\mathbb{C}\right\rbrace
在这上面计算形式幂级数 A(x) 的乘法逆元就不再需要 A(0)\neq 0 这个条件了，因为我们可以将 A(x) 写作 x^t\tilde{A}(x) ，我们先在形式幂级数环上计算 \tilde{A}(x) 的乘法逆元，最后再除以 x^t 就得到了 A(x) 在形式 Laurent 级数环上的乘法逆元了，这一步是自然的。
例 ：x^2 在形式 Laurent 级数环上的乘法逆元就是 x^{-2} 。
说点无关的东西，一些说法叫“减法卷积”的东西，其实就是形式 Laurent 级数的普通卷积而已（例如 Chirp Z 变换最早被称为 Fractional Fourier Transform 就是因为此）。
反（形式）Laurent 级数环
我们将 x 用 x^{-1} 替换就得到了 \mathbb{C}\left(\left(x^{-1}\right)\right) 。这么做的原因与我们在做多项式的带余除法时是一样的。在这个环上对于 Q(x)\in\mathbb{C}\left\lbrack x\right\rbrack 我们有
\begin{aligned}
\frac{1}{Q(x)}&=\cdots +f _ {0}x^{-\deg Q} \in \mathbb{C}\left(\left(x^{-1}\right)\right) \\
\frac{x^k}{Q(x)}&=\cdots +f _ kx^{-\deg Q}+\cdots +f _ 0x^{k-\deg Q} \in \mathbb{C}\left(\left(x^{-1}\right)\right) \\
\left\lbrack x^{\left(-\infty,0\right)}\right\rbrack \frac{x^k}{Q(x)}&=\left\lbrack x^{\left(-\infty,0\right)}\right\rbrack \frac{x^k\bmod{Q(x)}}{Q(x)}
\end{aligned}
而求出这些系数仍然可以使用翻转 Q(x) 再用形式幂级数的乘法逆元的算法来计算，因为
\begin{aligned}
\frac{1}{Q\left(x^{-1}\right)}=f_0x^{\deg Q}+\cdots\in\mathbb{C}\left(\left(x\right)\right) \\
\frac{1}{x^{\deg Q}Q\left(x^{-1}\right)}=f_0+\cdots\in\mathbb{C}\left\lbrack\left\lbrack x\right\rbrack\right\rbrack
\end{aligned}
也就是说所有这些翻转仅仅只是为了使得问题看起来简单，如果我们一定要在形式幂级数环 \mathbb{C}\left\lbrack\left\lbrack x\right\rbrack\right\rbrack 上计算所有结果也是可以的，只是表述起来比较麻烦。
Bostan–Mori 算法
我们这里只介绍原文中对应的 MSB-first 的算法来快速计算 x^k\bmod{Q(x)} 。根据前文我们已经明确实际上我们只需要计算
\left\lbrack x^{\left\lbrack -\deg Q,0\right)}\right\rbrack\frac{x^k}{Q(x)}
再和 Q(x) 作一次乘法，提取出 x^{\geq 0} 的系数就得到了 x^k\bmod{Q(x)} 。这是因为 x^{\lt -\deg Q} 的系数即便与 Q(x) 作乘法，其也只会贡献给 x^{\lt 0} 的项，不会影响到我们的结果。根据 Bostan–Mori 算法我们有
\frac{x^k}{Q(x)}=\frac{x^k}{Q(x)Q(-x)}\cdot Q(-x)
并且不论 k=2n 或 k=2n+1 ，我们都只需要计算
\left\lbrack x^{\left\lbrack -2\deg Q,0\right)}\right\rbrack \frac{x^{2n}}{Q(x)Q(-x)}
因为 Q(x)Q(-x) 是偶函数，所以 \left\lbrack x^{-2\deg Q-1}\right\rbrack \frac{x^{2n}}{Q(x)Q(-x)}=0 ，并且 x^{\lt -2\deg Q-1} 的系数与 Q(-x) 相乘之后只会贡献给 x^{\lt\deg Q-1} 的系数，在递归调用中是不会影响我们的返回值的，如果我们令 V\left(x^2\right)=Q(x)Q(-x) 那么我们只需计算
\left\lbrack x^{\left\lbrack -\deg Q,0\right)}\right\rbrack \frac{x^{n}}{V(x)}
当 k=0 时整个算法停止。
\begin{array}{ll}
&\textbf{Algorithm }\operatorname{\mathsf{RLSBostanMori}}\text{:} \\
&\textbf{Input}\text{: }Q(x),k\in\mathbb{N}\text{.} \\
&\textbf{Output}\text{: }\left\lbrack x^{\left\lbrack -\deg Q,0\right)}\right\rbrack \dfrac{x^k}{Q(x)},\text{where }Q(x)^{-1}\in\mathbb{C}\left(\left( x^{-1}\right)\right)\text{.} \\
1&\textbf{if }k=0 \textbf{ then return }\begin{bmatrix} \left(\left\lbrack x^{\deg Q}\right\rbrack Q(x)\right)^{-1} & 0 & \cdots & 0 \end{bmatrix} \\
2&V\left(x^2\right)\gets Q(x)Q(-x) \\
3&\begin{bmatrix} c _ {-\deg Q} & \cdots & c _ {-1}\end{bmatrix} \gets \operatorname{\mathsf{RLSBostanMori}}\left(V(x),\left\lfloor k/2\right\rfloor\right) \\
4&T(x)\gets \sum _ {j=0}^{-1+\deg Q}c _ {j-\deg Q}x^j \\
5&\sum _ {j=0}^{-1+3\deg Q} u _ jx^j\gets T\left(x^2\right)x^{k\bmod{2}}Q(-x) \\
6&\textbf{return }\begin{bmatrix} u _ {\deg Q} & \cdots & u _ {-1+2\deg Q}\end{bmatrix}
\end{array}
算法名中的 RLS 是反 Laurent 级数（Reversed Laurent Series）的英文缩写。
我们先对算法进行一些简单的说明，前三行都很简单，第四行我们为了让运算都是多项式，所以作了一定的处理，实际上可以直接考虑形如
\sum_{j=-\deg Q}^{-1}c_j x^{2j}
与 x^{k\bmod{2}}Q(-x) 作乘法，再提取乘积中的 x^{-\deg Q},\dots,x^{-1} 的系数。
后面我们再看另一种写法，当然我们的目标其实就是求出形式幂级数乘法逆元的一段连续系数，我当时用的一些更傻的方法来解释，现在可以给出一个几乎“对称”的解释了。比如给出形式 Laurent 级数 Q_0(x) ，我们希望求出 \left\lbrack x^{\left\lbrack k,k+\deg Q_0\right)}\right\rbrack Q_0(x)^{-1} ，也就是求出
\left\lbrack x^{\left\lbrack 0,\deg Q_0\right)}\right\rbrack \frac{x^{-k}}{Q_0(x)}
这个问题可以直接转换为上面的问题，但是直接处理呢？
\begin{array}{ll}
&\textbf{Algorithm }\operatorname{\mathsf{LSBostanMori}}\text{:} \\
&\textbf{Input}\text{: }Q_0(x),k\in\mathbb{N}\text{.} \\
&\textbf{Output}\text{: }\left\lbrack x^{\left\lbrack 0,\deg Q_0\right)}\right\rbrack \dfrac{x^{-k}}{Q_0(x)},\text{where }Q_0(x)^{-1}\in\mathbb{C}\left(\left( x\right)\right)\text{.} \\
1&\textbf{if }k=0 \textbf{ then return }\begin{bmatrix} \left\lbrack x^0\right\rbrack Q_0(x)^{-1} & \left\lbrack x^1\right\rbrack Q_0(x)^{-1} & \cdots & \left\lbrack x^{\deg Q_0-1}\right\rbrack Q_0(x)^{-1} \end{bmatrix} \\
2&V\left(x^2\right)\gets Q_0(x)Q_0(-x) \\
3&\begin{bmatrix} c _ {-\deg Q} & \cdots & c _ {-1}\end{bmatrix} \gets \operatorname{\mathsf{LSBostanMori}}\left(V(x),\left\lfloor k/2\right\rfloor\right) \\
4&T(x)\gets \sum _ {j=0}^{-1+\deg Q}c _ {j-\deg Q}x^j \\
5&\sum _ {j=0}^{-1+3\deg Q} u _ jx^j\gets T\left(x^2\right)x^{k\bmod{2}}Q(-x) \\
6&\textbf{return }\begin{bmatrix} u _ {\deg Q} & \cdots & u _ {-1+2\deg Q}\end{bmatrix}
\end{array}
还能优化吗？
当然可以，其实 Bostan–Mori 的论文已经给出了优化方法，但是我还是在这里随便说一下大概比较正确的调用方式。首先公认的事情是我们需要将上面给出伪代码的部分的输入输出的多项式都改成其 FFT 点值，这点毫无疑问，但是修改完毕之后这个算法仍然不能直接调用，而是需要在外层额外调用一次，然后对于其输出再进行单独的处理，细节也不是很多，可以看我的 代码和英文版的这篇文章大概。其中用到的 FFT 的一些处理技巧都在 FFT 里面描述了，大概都非常简单。
参考文献
Alin Bostan, Ryuhei Mori. A Simple and Fast Algorithm for Computing the N-th Term of a Linearly Recurrent Sequence. SOSA 2021: 118-132 url: https://arxiv.org/abs/2008.08822
Arne Storjohann. Algorithms for Matrix Canonical Forms. ETH Zürich. Diss., Technische Wissenschaften ETH Zürich, Nr. 13922, 2001. url: https://www.research-collection.ethz.ch/handle/20.500.11850/145127