跳转至

# 常系数齐次线性递推

## 简介

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

对于数列 (𝑎 𝑗)𝑗≥0(a j)j≥0 和其递推式

𝑎 𝑛=𝑑∑𝑗=1 𝑐 𝑗 𝑎 𝑛−𝑗,(𝑛≥𝑑)a n=∑j=1 d c j a n−j,(n≥d)

其中 𝑐 𝑗 c j 不全为零，我们的目标是在给出初值 𝑎 0,…,𝑎 𝑑−1 a 0,…,a d−1 和递推式中的 𝑐 1,…,𝑐 𝑑 c 1,…,c d 后求出 𝑎 𝑘 a k．如果 𝑘 ≫𝑑 k≫d，我们想要更快速的算法．

这里 (𝑎 𝑗)𝑗≥0(a j)j≥0 被称为 𝑑 d 阶的常系数齐次线性递推数列．

### Fiduccia 算法

Fiduccia 算法使用多项式取模和快速幂来计算 𝑎 𝑘 a k，时间为 𝑂(𝖬(𝑑)l o g⁡𝑘)O(M(d)log⁡k)，其中 𝑂(𝖬(𝑑))O(M(d)) 表示两个次数为 𝑂(𝑑)O(d) 的多项式相乘的时间．

算法：构造多项式 Γ(𝑥) :=𝑥 𝑑 −∑𝑑−1 𝑗=0 𝑐 𝑑−𝑗 𝑥 𝑗 Γ(x):=x d−∑j=0 d−1 c d−j x j 和 𝐴(𝑥) :=∑𝑑−1 𝑗=0 𝑎 𝑗 𝑥 𝑗 A(x):=∑j=0 d−1 a j x j，那么

𝑎 𝑘=⟨𝑥 𝑘 m o d Γ(𝑥),𝐴(𝑥)⟩a k=⟨x k mod Γ(x),A(x)⟩

其中定义 ⟨(∑𝑛−1 𝑗=0 𝑓 𝑗 𝑥 𝑗),(∑𝑛−1 𝑗=0 𝑔 𝑗 𝑥 𝑗)⟩ :=∑𝑛−1 𝑗=0 𝑓 𝑗 𝑔 𝑗⟨(∑j=0 n−1 f j x j),(∑j=0 n−1 g j x j)⟩:=∑j=0 n−1 f j g j 为内积．

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

𝐶 Γ:=⎡⎢ ⎢ ⎢ ⎢⎣𝑐 𝑑 1 𝑐 𝑑−1⋱⋮1 𝑐 1⎤⎥ ⎥ ⎥ ⎥⎦C Γ:=[c d 1 c d−1⋱⋮1 c 1]

我们定义多项式 𝑏(𝑥) :=∑𝑑−1 𝑗=0 𝑏 𝑗 𝑥 𝑗 b(x):=∑j=0 d−1 b j x j 和

𝐵 𝑏:=[𝑏 0 𝑏 1⋯𝑏 𝑑−1]⊺B b:=[b 0 b 1⋯b d−1]⊺

观察到

⎡⎢ ⎢ ⎢ ⎢⎣𝑐 𝑑 1 𝑐 𝑑−1⋱⋮1 𝑐 1⎤⎥ ⎥ ⎥ ⎥⎦⏟___⏟___⏟𝐶 Γ⎡⎢ ⎢ ⎢ ⎢⎣𝑏 0 𝑏 1⋮𝑏 𝑑−1⎤⎥ ⎥ ⎥ ⎥⎦⏟𝐵 𝑏=⎡⎢ ⎢ ⎢ ⎢⎣𝑐 𝑑 𝑏 𝑑−1 𝑏 0+𝑐 𝑑−1 𝑏 𝑑−1⋮𝑏 𝑑−2+𝑐 1 𝑏 𝑑−1⎤⎥ ⎥ ⎥ ⎥⎦⏟___⏟___⏟𝐵 𝑥 𝑏 m o d Γ[c d 1 c d−1⋱⋮1 c 1]⏟C Γ[b 0 b 1⋮b d−1]⏟B b=[c d b d−1 b 0+c d−1 b d−1⋮b d−2+c 1 b d−1]⏟B x b mod Γ

且

𝐶 Γ=[𝐵 𝑥 m o d Γ 𝐵 𝑥 2 m o d Γ⋯𝐵 𝑥 𝑑 m o d Γ],(𝐶 Γ)2=[𝐵 𝑥 2 m o d Γ 𝐵 𝑥 3 m o d Γ⋯𝐵 𝑥 𝑑+1 m o d Γ],⋯(𝐶 Γ)𝑘=[𝐵 𝑥 𝑘 m o d Γ 𝐵 𝑥 𝑘+1 m o d Γ⋯𝐵 𝑥 𝑘+𝑑 m o d Γ]C Γ=[B x mod Γ B x 2 mod Γ⋯B x d mod Γ],(C Γ)2=[B x 2 mod Γ B x 3 mod Γ⋯B x d+1 mod Γ],⋯(C Γ)k=[B x k mod Γ B x k+1 mod Γ⋯B x k+d mod Γ]

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

⎡⎢ ⎢ ⎢ ⎢⎣𝑎 𝑘 𝑎 𝑘+1⋮𝑎 𝑘+𝑑−1⎤⎥ ⎥ ⎥ ⎥⎦=⎡⎢ ⎢ ⎢ ⎢⎣1⋱1 𝑐 𝑑 𝑐 𝑑−1⋯𝑐 1⎤⎥ ⎥ ⎥ ⎥⎦𝑘⏟____⏟____⏟((𝐶 Γ)⊺)𝑘=((𝐶 Γ)𝑘)⊺⎡⎢ ⎢ ⎢ ⎢⎣𝑎 0 𝑎 1⋮𝑎 𝑑−1⎤⎥ ⎥ ⎥ ⎥⎦[a k a k+1⋮a k+d−1]=[1⋱1 c d c d−1⋯c 1]k⏟((C Γ)⊺)k=((C Γ)k)⊺[a 0 a 1⋮a d−1]

可知 ((𝐶 Γ)𝑘)⊺((C Γ)k)⊺ 的第一行为 𝐵 𝑥 𝑘 m o d Γ B x k mod Γ，根据矩阵乘法的定义得证．

### 表示为有理函数

对于上述数列 (𝑎 𝑗)𝑗≥0(a j)j≥0 一定存在有理函数

𝑃(𝑥)𝑄(𝑥)=∑𝑗≥0 𝑎 𝑗 𝑥 𝑗 P(x)Q(x)=∑j≥0 a j x j

且 𝑄(𝑥) =𝑥 𝑑 Γ(𝑥−1)Q(x)=x d Γ(x−1)，d e g⁡𝑃 <𝑑 deg⁡P<d．我们称其为「有理函数」是因为 𝑃(𝑥),𝑄(𝑥)P(x),Q(x) 是「多项式」．

证明：对于 𝑃(𝑥) =∑𝑑−1 𝑗=0 𝑝 𝑗 𝑥 𝑗 P(x)=∑j=0 d−1 p j x j 和 𝑄(𝑥) :=∑𝑑 𝑗=0 𝑞 𝑗 𝑥 𝑗 Q(x):=∑j=0 d q j x j 考虑 𝑃(𝑥)𝑄(𝑥) =∑𝑗≥0˜𝑞 𝑗 𝑥 𝑗 P(x)Q(x)=∑j≥0 q~j x j 的系数定义，这几乎就是形式幂级数「除法」的定义，

˜𝑞 𝑁=⎧{ { {⎨{ { {⎩𝑝 0 𝑞−1 0, i f 𝑁=0,(𝑝 𝑁−∑𝑁 𝑗=1 𝑞 𝑗˜𝑞 𝑁−𝑗)⋅𝑞−1 0, e l s e i f 𝑁<𝑑,−𝑞−1 0∑𝑑 𝑗=1 𝑞 𝑗˜𝑞 𝑁−𝑗, o t h e r w i s e.q~N={p 0 q 0−1, if N=0,(p N−∑j=1 N q j q~N−j)⋅q 0−1, else if N<d,−q 0−1∑j=1 d q j q~N−j, otherwise.

我们只需要令

𝑃(𝑥)=((∑𝑗≥0 𝑎 𝑗 𝑥 𝑗)⋅𝑥 𝑑 Γ(𝑥−1))m o d 𝑥 𝑑 P(x)=((∑j≥0 a j x j)⋅x d Γ(x−1))mod x d

那么根据 ˜𝑞 𝑁 q~N 的定义，必然有 𝑃(𝑥)𝑄(𝑥) =∑𝑗≥0 𝑎 𝑗 𝑥 𝑗 P(x)Q(x)=∑j≥0 a j x j．

### Bostan–Mori 算法

#### 计算单项

我们的目标仍然是给出上述多项式 𝑃(𝑥),𝑄(𝑥)P(x),Q(x)，求算 [𝑥 𝑘]𝑃(𝑥)𝑄(𝑥)[x k]P(x)Q(x)．

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

𝑃(𝑥)𝑄(𝑥)=𝑃(𝑥)𝑄(−𝑥)𝑄(𝑥)𝑄(−𝑥)=𝑈 0(𝑥 2)+𝑥 𝑈 1(𝑥 2)𝑉(𝑥 2)P(x)Q(x)=P(x)Q(−x)Q(x)Q(−x)=U 0(x 2)+x U 1(x 2)V(x 2)

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

[𝑥 𝑘]𝑃(𝑥)𝑄(𝑥)=[𝑥⌊𝑘/2⌋]𝑈 𝑘 m o d 2(𝑥)𝑉(𝑥)[x k]P(x)Q(x)=[x⌊k/2⌋]U k mod 2(x)V(x)

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

#### 计算连续若干项

目标是给出上述多项式 𝑃(𝑥),𝑄(𝑥)P(x),Q(x)，求算 [𝑥[𝐿,𝑅)]𝑃(𝑥)𝑄(𝑥)[x[L,R)]P(x)Q(x)．下面的计算中我们只需考虑对答案「有影响」的系数，这是 Bostan–Mori 算法的关键．

我们不妨假设 d e g⁡𝑃 <d e g⁡𝑄 deg⁡P<deg⁡Q，否则我们也可以通过一次带余除法使问题回到这种情况．

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

[𝑥[𝐿,𝑅)]1 𝑄(𝑥)=[𝑥[𝐿,𝑅)]1 𝑄(𝑥)𝑄(−𝑥)⋅𝑄(−𝑥)[x[L,R)]1 Q(x)=[x[L,R)]1 Q(x)Q(−x)⋅Q(−x)

我们需要求出 [𝑥[𝐿−d e g⁡𝑄,𝑅)]1 𝑄(𝑥)𝑄(−𝑥)[x[L−deg⁡Q,R)]1 Q(x)Q(−x) 然后作一次乘法并取出 𝑥 𝐿,…,𝑥 𝑅−1 x L,…,x R−1 的系数．令 𝑉(𝑥 2) =𝑄(𝑥)𝑄( −𝑥)V(x 2)=Q(x)Q(−x) 那么我们只需求出

[𝑥[⌈𝐿−d e g⁡𝑄 2⌉,⌈𝑅 2⌉)]1 𝑉(𝑥)[x[⌈L−deg⁡Q 2⌉,⌈R 2⌉)]1 V(x)

就可以还原出 [𝑥[𝐿−d e g⁡𝑄,𝑅)]1 𝑄(𝑥)𝑄(−𝑥)[x[L−deg⁡Q,R)]1 Q(x)Q(−x)．进而我们只需求出 [𝑥[𝐿−d e g⁡𝑃,𝑅)]1 𝑄(𝑥)[x[L−deg⁡P,R)]1 Q(x) 再和 𝑃(𝑥)P(x) 作一次乘法即可求出 [𝑥[𝐿,𝑅)]𝑃(𝑥)𝑄(𝑥)[x[L,R)]P(x)Q(x)．

上面的算法虽然已经可以工作，但是每一次的递归的时间复杂度与 𝑅 −𝐿 R−L 相关，我们希望能至少在递归求算时摆脱 𝑅 −𝐿 R−L，更具体的，我们先考虑求算 [𝑥[𝐿,𝐿+d e g⁡𝑄+1)]1 𝑄(𝑥)[x[L,L+deg⁡Q+1)]1 Q(x)，考虑

[𝑥[𝐿,𝐿+d e g⁡𝑄+1)]1 𝑄(𝑥)=[𝑥[𝐿,𝐿+d e g⁡𝑄+1)]1 𝑄(𝑥)𝑄(−𝑥)⋅𝑄(−𝑥)[x[L,L+deg⁡Q+1)]1 Q(x)=[x[L,L+deg⁡Q+1)]1 Q(x)Q(−x)⋅Q(−x)

我们需要求出

[𝑥[𝐿−d e g⁡𝑄,𝐿+d e g⁡𝑄+1)]1 𝑄(𝑥)𝑄(−𝑥)[x[L−deg⁡Q,L+deg⁡Q+1)]1 Q(x)Q(−x)

那么对于 𝑉(𝑥 2) =𝑄(𝑥)𝑄( −𝑥)V(x 2)=Q(x)Q(−x) 而言，我们只需求出

[𝑥[⌈(𝐿−d e g⁡𝑄)/2⌉,⌈(𝐿+d e g⁡𝑄+1)/2⌉)]1 𝑉(𝑥)[x[⌈(L−deg⁡Q)/2⌉,⌈(L+deg⁡Q+1)/2⌉)]1 V(x)

这是因为

[𝑥 𝑘]1 𝑄(𝑥)𝑄(−𝑥)=⎧{ {⎨{ {⎩[𝑥 𝑘/2]1 𝑉(𝑥),i f 𝑘≡0(m o d 2),0,o t h e r w i s e.[x k]1 Q(x)Q(−x)={[x k/2]1 V(x),if k≡0(mod 2),0,otherwise.

我们知道 𝐿 +d e g⁡𝑄 L+deg⁡Q 和 𝐿 −d e g⁡𝑄 L−deg⁡Q 的奇偶性是一样的，所以

⌈𝐿+d e g⁡𝑄+1 2⌉−⌈𝐿−d e g⁡𝑄 2⌉={d e g⁡𝑄+1,i f 𝐿+d e g⁡𝑄≡0(m o d 2),d e g⁡𝑄,o t h e r w i s e.⌈L+deg⁡Q+1 2⌉−⌈L−deg⁡Q 2⌉={deg⁡Q+1,if L+deg⁡Q≡0(mod 2),deg⁡Q,otherwise.

这样我们可以写出伪代码

𝐀 𝐥 𝐠 𝐨 𝐫 𝐢 𝐭 𝐡 𝐦 S l i c e-C o e f f i c i e n t s⁡(𝑄,𝐿):𝐈 𝐧 𝐩 𝐮 𝐭: 𝑄(𝑥)∈ℂ[𝑥],𝐿∈ℤ.𝐎 𝐮 𝐭 𝐩 𝐮 𝐭: [𝑥[𝐿,𝐿+d e g⁡𝑄+1)]𝑄(𝑥)−1.1 𝐢 𝐟 𝐿≤1 𝐭 𝐡 𝐞 𝐧 𝐫 𝐞 𝐭 𝐮 𝐫 𝐧 [𝑥[𝐿,𝐿+d e g⁡𝑄+1)]𝑄(𝑥)−1 U s e o t h e r a l g o r i t h m t o c o m p u t e 𝑄(𝑥)−1 2 𝑉(𝑥 2)←𝑄(𝑥)𝑄(−𝑥)3 𝑘←⌈𝐿−d e g⁡𝑄 2⌉4(𝑡 𝑘,…,𝑡 𝑘+d e g⁡𝑄)←S l i c e-C o e f f i c i e n t s⁡(𝑉,𝑘)5 𝑇(𝑥)←𝑥(𝐿−d e g⁡𝑄)m o d 2∑d e g⁡𝑄 𝑗=0 𝑡 𝑗+𝑘 𝑥 2 𝑗 6 𝐫 𝐞 𝐭 𝐮 𝐫 𝐧 [𝑥[d e g⁡𝑄,2 d e g⁡𝑄+1)]𝑇(𝑥)𝑄(−𝑥)Algorithm Slice-Coefficients⁡(Q,L):Input: Q(x)∈C[x],L∈Z.Output: [x[L,L+deg⁡Q+1)]Q(x)−1.1 if L≤1 then return [x[L,L+deg⁡Q+1)]Q(x)−1 Use other algorithm to compute Q(x)−1 2 V(x 2)←Q(x)Q(−x)3 k←⌈L−deg⁡Q 2⌉4(t k,…,t k+deg⁡Q)←Slice-Coefficients⁡(V,k)5 T(x)←x(L−deg⁡Q)mod 2∑j=0 deg⁡Q t j+k x 2 j 6 return [x[deg⁡Q,2 deg⁡Q+1)]T(x)Q(−x)

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

#### 找到新的有理函数表示

我们知道 𝑄(𝑥)Q(x) 本身和 𝑄(𝑥)−1 Q(x)−1 的一部分连续的系数比如 [𝑥[𝐿,𝐿+d e g⁡𝑄)]𝑄(𝑥)−1[x[L,L+deg⁡Q)]Q(x)−1 和 𝐿 ≥0 L≥0，我们希望求出 [𝑥[𝐿+d e g⁡𝑄,𝐿+2 d e g⁡𝑄)]𝑄(𝑥)−1[x[L+deg⁡Q,L+2 deg⁡Q)]Q(x)−1，这等价于我们要求某个 𝑃(𝑥)P(x) 且 d e g⁡𝑃 <d e g⁡𝑄 deg⁡P<deg⁡Q 使得 𝑃(𝑥)𝑄(𝑥)P(x)Q(x) 的前 d e g⁡𝑄 deg⁡Q 项与 [𝑥[𝐿,𝐿+d e g⁡𝑄)]𝑄(𝑥)−1[x[L,L+deg⁡Q)]Q(x)−1 相同．简单来说：递推关系（有理函数的分母）是不变的，我们所做的只是更换初值（有理函数的分子）．

具体的，考虑

𝑃(𝑥)𝑄(𝑥)=∑𝑗≥0 𝑎 𝑗 𝑥 𝑗 P(x)Q(x)=∑j≥0 a j x j

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

∑𝑗≥𝑛 𝑎 𝑗 𝑥 𝑗−𝑛=𝑃(𝑥)𝑄(𝑥)𝑥 𝑛−𝑄(𝑥)∑𝑛−1 𝑗=0 𝑎 𝑗 𝑥 𝑗 𝑄(𝑥)𝑥 𝑛∑j≥n a j x j−n=P(x)Q(x)x n−Q(x)∑j=0 n−1 a j x j Q(x)x n

我们先用一次 S l i c e-C o e f f i c i e n t s⁡(𝑄,𝐿 −d e g⁡𝑃)Slice-Coefficients⁡(Q,L−deg⁡P) 计算出 [𝑥[𝐿−d e g⁡𝑃,𝐿−d e g⁡𝑃+d e g⁡𝑄+1)]𝑄(𝑥)−1[x[L−deg⁡P,L−deg⁡P+deg⁡Q+1)]Q(x)−1，然后我们扩展合并出 [𝑥[𝐿−d e g⁡𝑃,𝐿+d e g⁡𝑄)]𝑄(𝑥)−1[x[L−deg⁡P,L+deg⁡Q)]Q(x)−1，再重新计算一个分子使得

̃𝑃(𝑥)𝑄(𝑥)=∑𝑗≥0([𝑥 𝐿+𝑗]𝑃(𝑥)𝑄(𝑥))𝑥 𝑗 P~(x)Q(x)=∑j≥0([x L+j]P(x)Q(x))x j

最后我们使用形式幂级数的除法计算出 [𝑥[0,𝑅−𝐿)]̃𝑃(𝑥)𝑄(𝑥)[x[0,R−L)]P~(x)Q(x)，时间为 𝑂(𝖬(𝑑)l o g⁡𝐿 +𝖬(𝑅 −𝐿))O(M(d)log⁡L+M(R−L))．

## 参考文献

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

> 本页面最近更新：2026/1/7 08:56:54，更新历史
> 发现错误？想一起完善？ 在 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, qz-cqy, SamZhangQingChuan, sshwy, StudyingFather, test12345-pupil, thredreams, TrisolarisHD, untitledunrevised, Xeonacid
> 本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用