常系数齐次线性递推

简介

常系数齐次线性递推数列（又称为 C-finite 或 C-recursive 数列）是常见的一类基础的递推数列。

对于数列 ${\left(a_j\right)}_{j\ge 0}$ 和其递推式

$$a_n=\sum _{j=1}^d{c}_j{a}_{n-j},(n\ge d)$$

其中 $c_j$ 不全为零，我们的目标是在给出初值 $a_0,\dots ,a_{d-1}$ 和递推式中的 $c_1,\dots ,c_d$ 后求出 $a_k$。如果 $k\gg d$，我们想要更快速的算法。

这里 ${\left(a_j\right)}_{j\ge 0}$ 被称为 $d$ 阶的常系数齐次线性递推数列。

Fiduccia 算法

Fiduccia 算法使用多项式取模和快速幂来计算 $a_k$，时间为 $O(\mathsf{M}(d)\log k)$，其中 $O(\mathsf{M}(d))$ 表示两个次数为 $O(d)$ 的多项式相乘的时间。

算法：构造多项式 $\mathrm{\Gamma }(x):=x^d-\sum _{j=0}^{d-1}{c}_{d-j}{x}^j$ 和 $A(x):=\sum _{j=0}^{d-1}{a}_j{x}^j$，那么

$$a_k=⟨x^{k}mod\mathrm{\Gamma }(x),A(x)⟩$$

其中定义 $⟨\left(\sum _{j=0}^{n-1}{f}_j{x}^j\right),\left(\sum _{j=0}^{n-1}{g}_j{x}^j\right)⟩:=\sum _{j=0}^{n-1}{f}_j{g}_j$ 为内积。

证明：我们定义 $\mathrm{\Gamma }(x)$ 的友矩阵为

$$C_{\mathrm{\Gamma }}:=\left[\begin{array}{cccc}& & & c_d\\ 1& & & c_{d-1}\\ & \ddots & & ⋮\\ & & 1& c_1\end{array}\right]$$

我们定义多项式 $b(x):=\sum _{j=0}^{d-1}{b}_j{x}^j$ 和

$$B_b:={\left[\begin{array}{cccc}{b}_0& b_1& \cdots & b_{d-1}\end{array}\right]}^{⊺}$$

观察到

$$\underset{C_{\mathrm{\Gamma }}}{\underbrace{\left[\begin{array}{cccc}& & & c_d\\ 1& & & c_{d-1}\\ & \ddots & & ⋮\\ & & 1& c_1\end{array}\right]}}\underset{B_b}{\underbrace{\left[\begin{array}{c}{b}_0\\ b_1\\ ⋮\\ b_{d-1}\end{array}\right]}}=\underset{B_{xbmod\mathrm{\Gamma }}}{\underbrace{\left[\begin{array}{c}{c}_d{b}_{d-1}\\ b_0+c_{d-1}{b}_{d-1}\\ ⋮\\ b_{d-2}+c_1{b}_{d-1}\end{array}\right]}}$$

且

$$\begin{array}{rl}{C}_{\mathrm{\Gamma }}& =\left[\begin{array}{cccc}{B}_{xmod\mathrm{\Gamma }}& B_{x^{2}mod\mathrm{\Gamma }}& \cdots & B_{x^{d}mod\mathrm{\Gamma }}\end{array}\right],\\ {\left(C_{\mathrm{\Gamma }}\right)}^2& =\left[\begin{array}{cccc}{B}_{x^{2}mod\mathrm{\Gamma }}& B_{x^{3}mod\mathrm{\Gamma }}& \cdots & B_{x^{d+1}mod\mathrm{\Gamma }}\end{array}\right],\\ \cdots \\ {\left(C_{\mathrm{\Gamma }}\right)}^k& =\left[\begin{array}{cccc}{B}_{x^{k}mod\mathrm{\Gamma }}& B_{x^{k+1}mod\mathrm{\Gamma }}& \cdots & B_{x^{k+d}mod\mathrm{\Gamma }}\end{array}\right]\end{array}$$

我们将这个递推用矩阵表示有

$$\left[\begin{array}{c}{a}_k\\ a_{k+1}\\ ⋮\\ a_{k+d-1}\end{array}\right]=\underset{{\left({\left(C_{\mathrm{\Gamma }}\right)}^{⊺}\right)}^k={\left({\left(C_{\mathrm{\Gamma }}\right)}^k\right)}^{⊺}}{\underbrace{{\left[\begin{array}{cccc}& 1& & \\ & & \ddots & \\ & & & 1\\ c_d& c_{d-1}& \cdots & c_1\end{array}\right]}^k}}\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_{d-1}\end{array}\right]$$

可知 ${\left({\left(C_{\mathrm{\Gamma }}\right)}^k\right)}^{⊺}$ 的第一行为 $B_{x^{k}mod\mathrm{\Gamma }}$，根据矩阵乘法的定义得证。

表示为有理函数

对于上述数列 ${\left(a_j\right)}_{j\ge 0}$ 一定存在有理函数

$$\frac{P(x)}{Q(x)}=\sum _{j\ge 0}{a}_j{x}^j$$

且 $Q(x)=x^d\mathrm{\Gamma }\left(x^{-1}\right)$，$\mathrm{deg}P<d$。我们称其为「有理函数」是因为 $P(x),Q(x)$ 是「多项式」。

证明：对于 $P(x)=\sum _{j=0}^{d-1}{p}_j{x}^j$ 和 $Q(x):=\sum _{j=0}^d{q}_j{x}^j$ 考虑 ${\displaystyle \frac{P(x)}{Q(x)}}=\sum _{j\ge 0}{\tilde{q}}_j{x}^j$ 的系数定义，这几乎就是形式幂级数「除法」的定义，

$${\tilde{q}}_N=\{\begin{array}{ll}{p}_0{q}_0^{-1},& \text{ if }N=0,\\ (p_N-\sum _{j=1}^N{q}_j{\tilde{q}}_{N-j})\cdot q_0^{-1},& \text{ else if }N<d,\\ -q_0^{-1}\sum _{j=1}^d{q}_j{\tilde{q}}_{N-j},& \text{ otherwise}.\end{array}$$

我们只需要令

$$P(x)=(\left(\sum _{j\ge 0}{a}_j{x}^j\right)\cdot x^d\mathrm{\Gamma }\left(x^{-1}\right))mod{x}^d$$

那么根据 ${\tilde{q}}_N$ 的定义，必然有 ${\displaystyle \frac{P(x)}{Q(x)}}=\sum _{j\ge 0}{a}_j{x}^j$。

Bostan–Mori 算法

计算单项

我们的目标仍然是给出上述多项式 $P(x),Q(x)$，求算 $\left[x^k\right]{\displaystyle \frac{P(x)}{Q(x)}}$。

Bostan–Mori 算法基于 Graeffe 迭代，对于上述多项式 $P(x),Q(x)$ 有

$$\frac{P(x)}{Q(x)}=\frac{P(x)Q(-x)}{Q(x)Q(-x)}=\frac{U_0(x^2)+x{U}_1(x^2)}{V(x^2)}$$

因为分母 $V(x^2)$ 是偶函数，所以子问题只需考虑其中的一侧

$$\left[x^k\right]{\displaystyle \frac{P(x)}{Q(x)}}=\left[x^{\lfloor k/2\rfloor }\right]\frac{U_{kmod2}(x)}{V(x)}$$

我们付出两次多项式乘法的代价使得问题至少减少为原先的一半，而当 $k=0$ 时显然有 $\left[x^0\right]{\displaystyle \frac{P(x)}{Q(x)}}={\displaystyle \frac{P(0)}{Q(0)}}$，时间复杂度同上。

计算连续若干项

目标是给出上述多项式 $P(x),Q(x)$，求算 $\left[x^{[L,R)}\right]{\displaystyle \frac{P(x)}{Q(x)}}$。下面的计算中我们只需考虑对答案「有影响」的系数，这是 Bostan–Mori 算法的关键。

我们不妨假设 $\mathrm{deg}P<\mathrm{deg}Q$，否则我们也可以通过一次带余除法使问题回到这种情况。

我们先考虑更简单的问题：

$$\left[x^{[L,R)}\right]\frac{1}{Q(x)}=\left[x^{[L,R)}\right]\frac{1}{Q(x)Q(-x)}\cdot Q(-x)$$

我们需要求出 $\left[x^{[L-\mathrm{deg}Q,R)}\right]{\displaystyle \frac{1}{Q(x)Q(-x)}}$ 然后作一次乘法并取出 $x^L,\dots ,x^{R-1}$ 的系数。令 $V(x^2)=Q(x)Q(-x)$ 那么我们只需求出

$$\left[x^{[\lceil \frac{L-\mathrm{deg}Q}{2}\rceil ,\lceil \frac{R}{2}\rceil )}\right]\frac{1}{V(x)}$$

就可以还原出 $\left[x^{[L-\mathrm{deg}Q,R)}\right]{\displaystyle \frac{1}{Q(x)Q(-x)}}$。进而我们只需求出 $\left[x^{[L-\mathrm{deg}P,R)}\right]{\displaystyle \frac{1}{Q(x)}}$ 再和 $P(x)$ 作一次乘法即可求出 $\left[x^{[L,R)}\right]{\displaystyle \frac{P(x)}{Q(x)}}$。

上面的算法虽然已经可以工作，但是每一次的递归的时间复杂度与 $R-L$ 相关，我们希望能至少在递归求算时摆脱 $R-L$，更具体的，我们先考虑求算 $\left[x^{[L,L+\mathrm{deg}Q+1)}\right]{\displaystyle \frac{1}{Q(x)}}$，考虑

$$\left[x^{[L,L+\mathrm{deg}Q+1)}\right]\frac{1}{Q(x)}=\left[x^{[L,L+\mathrm{deg}Q+1)}\right]{\displaystyle \frac{1}{Q(x)Q(-x)}}\cdot Q(-x)$$

我们需要求出

$$\left[x^{[L-\mathrm{deg}Q,L+\mathrm{deg}Q+1)}\right]{\displaystyle \frac{1}{Q(x)Q(-x)}}$$

那么对于 $V(x^2)=Q(x)Q(-x)$ 而言，我们只需求出

$$\left[x^{[\lceil (L-\mathrm{deg}Q)/2\rceil ,\lceil (L+\mathrm{deg}Q+1)/2\rceil )}\right]\frac{1}{V(x)}$$

这是因为

$$\left[x^k\right]{\displaystyle \frac{1}{Q(x)Q(-x)}}=\{\begin{array}{ll}\left[x^{k/2}\right]{\displaystyle \frac{1}{V(x)}},& \text{if }k\equiv 0(\mathrm{mod}2),\\ 0,& \text{otherwise}.\end{array}$$

我们知道 $L+\mathrm{deg}Q$ 和 $L-\mathrm{deg}Q$ 的奇偶性是一样的，所以

$$\lceil \frac{L+\mathrm{deg}Q+1}{2}\rceil -\lceil \frac{L-\mathrm{deg}Q}{2}\rceil =\{\begin{array}{ll}\mathrm{deg}Q+1,& \text{if }L+\mathrm{deg}Q\equiv 0(\mathrm{mod}2),\\ \mathrm{deg}Q,& \text{otherwise}.\end{array}$$

这样我们可以写出伪代码

$$\begin{array}{ll}& \mathbf{\text{Algorithm }}\mathrm{Slice-Coefficients}(Q,L)\text{:}\\ & \mathbf{\text{Input}}\text{: }Q(x)\in \mathbb{C}\left[x\right],L\in \mathbb{Z}\text{.}\\ & \mathbf{\text{Output}}\text{: }\left[x^{[L,L+\mathrm{deg}Q+1)}\right]Q(x{)}^{-1}\text{.}\\ 1& \mathbf{\text{if }}L\le 1\mathbf{\text{ then return }}\left[x^{[L,L+\mathrm{deg}Q+1)}\right]Q(x{)}^{-1}\\ & \text{Use other algorithm to compute }Q(x{)}^{-1}\\ 2& V(x^2)\leftarrow Q(x)Q(-x)\\ 3& k\leftarrow \lceil \frac{L-\mathrm{deg}Q}{2}\rceil \\ 4& (t_k,\dots ,t_{k+\mathrm{deg}Q})\leftarrow \mathrm{Slice-Coefficients}(V,k)\\ 5& T(x)\leftarrow x^{(L-\mathrm{deg}Q)mod2}\sum _{j=0}^{\mathrm{deg}Q}{t}_{j+k}{x}^{2j}\\ 6& \mathbf{\text{return }}\left[x^{[\mathrm{deg}Q,2\mathrm{deg}Q+1)}\right]T(x)Q(-x)\end{array}$$

但是只有这个算法还不够，我们需要重新找到一个有理函数并求算更多系数。

找到新的有理函数表示

我们知道 $Q(x)$ 本身和 $Q(x{)}^{-1}$ 的一部分连续的系数比如 $\left[x^{[L,L+\mathrm{deg}Q)}\right]Q(x{)}^{-1}$ 和 $L\ge 0$，我们希望求出 $\left[x^{[L+\mathrm{deg}Q,L+2\mathrm{deg}Q)}\right]Q(x{)}^{-1}$，这等价于我们要求某个 $P(x)$ 且 $\mathrm{deg}P<\mathrm{deg}Q$ 使得 ${\displaystyle \frac{P(x)}{Q(x)}}$ 的前 $\mathrm{deg}Q$ 项与 $\left[x^{[L,L+\mathrm{deg}Q)}\right]Q(x{)}^{-1}$ 相同。简单来说：递推关系（有理函数的分母）是不变的，我们所做的只是更换初值（有理函数的分子）。

具体的，考虑

$$\frac{P(x)}{Q(x)}=\sum _{j\ge 0}{a}_j{x}^j$$

我们现在希望将递推前进 $n$ 项，那么就是

$$\sum _{j\ge n}{a}_j{x}^{j-n}=\frac{P(x)}{Q(x)x^n}-\frac{Q(x)\sum _{j=0}^{n-1}{a}_j{x}^j}{Q(x)x^n}$$

我们先用一次 $\mathrm{Slice-Coefficients}(Q,L-\mathrm{deg}P)$ 计算出 $\left[x^{[L-\mathrm{deg}P,L-\mathrm{deg}P+\mathrm{deg}Q+1)}\right]Q(x{)}^{-1}$，然后我们扩展合并出 $\left[x^{[L-\mathrm{deg}P,L+\mathrm{deg}Q)}\right]Q(x{)}^{-1}$，再重新计算一个分子使得

$$\frac{\tilde{P}(x)}{Q(x)}=\sum _{j\ge 0}\left(\left[x^{L+j}\right]\frac{P(x)}{Q(x)}\right)x^j$$

最后我们使用形式幂级数的除法计算出 $\left[x^{[0,R-L)}\right]{\displaystyle \frac{\tilde{P}(x)}{Q(x)}}$，时间为 $O(\mathsf{M}(d)\log L+\mathsf{M}(R-L))$。

参考文献

1. Alin Bostan, Ryuhei Mori.A Simple and Fast Algorithm for Computing the $N$-th Term of a Linearly Recurrent Sequence.

本页面最近更新：2024/10/30 21:50:12，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：Tiphereth-A, Enter-tainer, hly1204, QAQAutoMaton, Marcythm, 97littleleaf11, ksyx, shuzhouliu, abc1763613206, c-forrest, CCXXXI, Early0v0, EndlessCheng, fps5283, Great-designer, H-J-Granger, hsfzLZH1, huayucaiji, Ir1d, kenlig, ouuan, ouuan, SamZhangQingChuan, sshwy, StudyingFather, test12345-pupil, thredreams, TrisolarisHD, untitledunrevised, Xeonacid
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用