常系数齐次线性递推

简介

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

对于数列 ${(a_{j})}_{j \geq 0}$ 和其递推式

a_{n} = \sum_{j = 1}^{d} c_{j} a_{n - j}, (n \geq d)

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

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

Fiduccia 算法

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

算法：构造多项式 $Γ (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 Γ (x), A (x) ⟩

其中定义 $⟨ (\sum_{j = 0}^{n - 1} f_{j} x^{j}), (\sum_{j = 0}^{n - 1} g_{j} x^{j}) ⟩ := \sum_{j = 0}^{n - 1} f_{j} g_{j}$ 为内积。

证明：我们定义 $Γ (x)$ 的友矩阵为

C_{Γ} := [\begin{matrix} c_{d} \\ 1 & c_{d - 1} \\ ⋱ & ⋮ \\ 1 & c_{1} \end{matrix}]

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

B_{b} := {[\begin{matrix} b_{0} & b_{1} & \dots & b_{d - 1} \end{matrix}]}^{⊺}

观察到

\underset{C_{Γ}}{\underset{⏟}{[\begin{matrix} c_{d} \\ 1 & c_{d - 1} \\ ⋱ & ⋮ \\ 1 & c_{1} \end{matrix}]}} \underset{B_{b}}{\underset{⏟}{[\begin{matrix} b_{0} \\ b_{1} \\ ⋮ \\ b_{d - 1} \end{matrix}]}} = \underset{B_{x b mod Γ}}{\underset{⏟}{[\begin{matrix} c_{d} b_{d - 1} \\ b_{0} + c_{d - 1} b_{d - 1} \\ ⋮ \\ b_{d - 2} + c_{1} b_{d - 1} \end{matrix}]}}

且

\begin{aligned} C_{Γ} & = [\begin{matrix} B_{x mod Γ} & B_{x^{2} mod Γ} & \dots & B_{x^{d} mod Γ} \end{matrix}], \\ {(C_{Γ})}^{2} & = [\begin{matrix} B_{x^{2} mod Γ} & B_{x^{3} mod Γ} & \dots & B_{x^{d + 1} mod Γ} \end{matrix}], \\ \dots \\ {(C_{Γ})}^{k} & = [\begin{matrix} B_{x^{k} mod Γ} & B_{x^{k + 1} mod Γ} & \dots & B_{x^{k + d} mod Γ} \end{matrix}] \end{aligned}

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

[\begin{matrix} a_{k} \\ a_{k + 1} \\ ⋮ \\ a_{k + d - 1} \end{matrix}] = \underset{{({(C_{Γ})}^{⊺})}^{k} = {({(C_{Γ})}^{k})}^{⊺}}{\underset{⏟}{{[\begin{matrix} 1 \\ ⋱ \\ 1 \\ c_{d} & c_{d - 1} & \dots & c_{1} \end{matrix}]}^{k}}} [\begin{matrix} a_{0} \\ a_{1} \\ ⋮ \\ a_{d - 1} \end{matrix}]

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

表示为有理函数

对于上述数列 ${(a_{j})}_{j \geq 0}$ 一定存在有理函数

\frac{P (x)}{Q (x)} = \sum_{j \geq 0} a_{j} x^{j}

且 $Q (x) = x^{d} Γ (x^{- 1})$ ， $\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}$ 考虑 $\frac{P (x)}{Q (x)} = \sum_{j \geq 0} {\tilde{q}}_{j} x^{j}$ 的系数定义，这几乎就是形式幂级数「除法」的定义，

{\tilde{q}}_{N} = {\begin{cases} p_{0} q_{0}^{- 1}, & if N = 0, \\ (p_{N} - \sum_{j = 1}^{N} q_{j} {\tilde{q}}_{N - j}) \cdot q_{0}^{- 1}, & else if N < d, \\ - q_{0}^{- 1} \sum_{j = 1}^{d} q_{j} {\tilde{q}}_{N - j}, & otherwise . \end{cases}

我们只需要令

P (x) = ((\sum_{j \geq 0} a_{j} x^{j}) \cdot x^{d} Γ (x^{- 1})) mod x^{d}

那么根据 ${\tilde{q}}_{N}$ 的定义，必然有 $\frac{P (x)}{Q (x)} = \sum_{j \geq 0} a_{j} x^{j}$ 。

Bostan–Mori 算法

计算单项

我们的目标仍然是给出上述多项式 $P (x), Q (x)$ ，求算 $[x^{k}] \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})$ 是偶函数，所以子问题只需考虑其中的一侧

[x^{k}] \frac{P (x)}{Q (x)} = [x^{⌊ k / 2 ⌋}] \frac{U_{k mod 2} (x)}{V (x)}

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

计算连续若干项

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

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

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

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

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

[x^{[⌈ \frac{L - \deg Q}{2} ⌉, ⌈ \frac{R}{2} ⌉)}] \frac{1}{V (x)}

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

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

[x^{[L, L + \deg Q + 1)}] \frac{1}{Q (x)} = [x^{[L, L + \deg Q + 1)}] \frac{1}{Q (x) Q (- x)} \cdot Q (- x)

我们需要求出

[x^{[L - \deg Q, L + \deg Q + 1)}] \frac{1}{Q (x) Q (- x)}

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

[x^{[⌈ (L - \deg Q) / 2 ⌉, ⌈ (L + \deg Q + 1) / 2 ⌉)}] \frac{1}{V (x)}

这是因为

[x^{k}] \frac{1}{Q (x) Q (- x)} = {\begin{cases} [x^{k / 2}] \frac{1}{V (x)}, & if k \equiv 0 (\mod 2), \\ 0, & otherwise . \end{cases}

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

⌈ \frac{L + \deg Q + 1}{2} ⌉ - ⌈ \frac{L - \deg Q}{2} ⌉ = {\begin{cases} \deg Q + 1, & if L + \deg Q \equiv 0 (\mod 2), \\ \deg Q, & otherwise . \end{cases}

这样我们可以写出伪代码

\begin{array}{ll} Algorithm Slice-Coefficients (Q, L) : \\ Input : Q (x) \in C [x], L \in Z . \\ Output : [x^{[L, L + \deg Q + 1)}] Q (x)^{- 1} . \\ 1 & if L \leq 1 then return [x^{[L, L + \deg Q + 1)}] Q (x)^{- 1} \\ Use other algorithm to compute Q (x)^{- 1} \\ 2 & V (x^{2}) \leftarrow Q (x) Q (- x) \\ 3 & k \leftarrow ⌈ \frac{L - \deg Q}{2} ⌉ \\ 4 & (t_{k}, \dots, t_{k + \deg Q}) \leftarrow Slice-Coefficients (V, k) \\ 5 & T (x) \leftarrow x^{(L - \deg Q) mod 2} \sum_{j = 0}^{\deg Q} t_{j + k} x^{2 j} \\ 6 & return [x^{[\deg Q, 2 \deg Q + 1)}] T (x) Q (- x) \end{array}

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

找到新的有理函数表示

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

具体的，考虑

\frac{P (x)}{Q (x)} = \sum_{j \geq 0} a_{j} x^{j}

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

\sum_{j \geq 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}}

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

\frac{\tilde{P} (x)}{Q (x)} = \sum_{j \geq 0} ([x^{L + j}] \frac{P (x)}{Q (x)}) x^{j}

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

参考文献

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 协议之条款下提供，附加条款亦可能应用