Linear Recurrence with Constant Coefficients

Introduction

A constant-coefficient homogeneous linear recurrence sequence (also known as C-finite or C-recursive sequence) is a common type of fundamental recurrence sequence.

For a sequence \(\left(a_j\right)_{j\geq 0}\) and its recurrence relation

\[ a_n=\sum_{j=1}^{d}c_ja_{n-j},\qquad (n\geq d) \]

where not all \(c_j\) are zero, our goal is to compute \(a_k\) given the initial values \(a_0,\dots ,a_{d-1}\) and the coefficients \(c_1,\dots ,c_d\) in the recurrence. If \(k\gg d\), we desire a faster algorithm.

Here, \(\left(a_j\right)_{j\geq 0}\) is called a constant-coefficient homogeneous linear recurrence sequence of order \(d\).

Fiduccia's Algorithm

Fiduccia's algorithm uses polynomial modular reduction and fast exponentiation to compute \(a_k\) in \(O(\mathsf{M}(d)\log k)\) time, where \(O(\mathsf{M}(d))\) denotes the time to multiply two polynomials of degree \(O(d)\).

Algorithm: Construct the polynomials \(\Gamma(x):=x^d-\sum_{j=0}^{d-1}c_{d-j}x^j\) and \(A(x):=\sum_{j=0}^{d-1}a_jx^j\). Then

\[ a_k=\left\langle x^k\bmod{\Gamma(x)},A(x)\right\rangle \]

where the inner product is defined as \(\left\langle \left(\sum_{j=0}^{n-1}f_jx^j\right),\left(\sum_{j=0}^{n-1}g_jx^j\right) \right\rangle :=\sum_{j=0}^{n-1}f_jg_j\).

Proof: We define the companion matrix of \(\Gamma(x)\) as

\[ C_\Gamma:= \begin{bmatrix} &&&c_d \\ 1&&&c_{d-1} \\ &\ddots &&\vdots \\ &&1&c_1 \end{bmatrix} \]

We define a polynomial \(b(x):=\sum_{j=0}^{d-1}b_jx^j\) and

\[ B_b:=\begin{bmatrix}b_0&b_1&\cdots &b_{d-1}\end{bmatrix}^{\intercal} \]

Observe that

\[ \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_db_{d-1} \\ b_0+c_{d-1}b_{d-1} \\ \vdots \\ b_{d-2}+c_1b_{d-1} \end{bmatrix}} _ {B_{xb\bmod{\Gamma}}} \]

and

\[ \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}, \\ \cdots \\ \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} \]

Representing the recurrence in matrix form gives

\[ \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} \]

It follows that the first row of \(\left(\left(C_\Gamma\right)^{k}\right)^{\intercal}\) is \(B_{x^k\bmod{\Gamma}}\), proving the claim by the definition of matrix multiplication.

Representation as a Rational Function

For the above sequence \(\left(a_j\right)_{j\geq 0}\), there always exists a rational function

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

with \(Q(x)=x^d\Gamma\left(x^{-1}\right)\) and \(\deg{P}<d\). We call it a "rational function" because \(P(x)\) and \(Q(x)\) are "polynomials".

Proof: For \(P(x)=\sum_{j=0}^{d-1}p_jx^j\) and \(Q(x):=\sum_{j=0}^{d}q_jx^j\), consider the coefficient definition of \(\dfrac{P(x)}{Q(x)}=\sum_{j\geq 0}\tilde{q}_jx^j\). This is essentially the definition of "division" for formal power series:

\[ \tilde{q}_N= \begin{cases} p_0q_0^{-1},&\text{ if }N=0, \\ \left(p_N-\sum_{j=1}^{N}q_j\tilde{q}_{N-j}\right)\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{cases} \]

We only need to set

\[ P(x)=\left(\left(\sum_{j\geq 0}a_jx^j\right)\cdot x^d\Gamma\left(x^{-1}\right)\right)\bmod{x^d} \]

Then, by the definition of \(\tilde{q}_N\), we must have \(\dfrac{P(x)}{Q(x)}=\sum_{j\geq 0}a_jx^j\).

Bostan–Mori Algorithm

Computing a Single Term

Our goal is still, given the polynomials \(P(x),Q(x)\) above, to compute \(\left\lbrack x^k\right\rbrack\dfrac{P(x)}{Q(x)}\).

The Bostan–Mori algorithm is based on Graeffe iteration. For the polynomials \(P(x),Q(x)\) above, we have

\[ \frac{P(x)}{Q(x)}=\frac{P(x)Q(-x)}{Q(x)Q(-x)}=\frac{U_0(x^2)+xU_1(x^2)}{V(x^2)} \]

Since the denominator \(V(x^2)\) is an even function, the subproblem only needs to consider one side:

\[ \left\lbrack x^k\right\rbrack\dfrac{P(x)}{Q(x)}=\left\lbrack x^{\left\lfloor k/2\right\rfloor}\right\rbrack \frac{U_{k\bmod{2}}(x)}{V(x)} \]

We pay the cost of two polynomial multiplications to reduce the problem size by at least half. When \(k=0\), clearly \(\left\lbrack x^0\right\rbrack \dfrac{P(x)}{Q(x)}=\dfrac{P(0)}{Q(0)}\). The time complexity is the same as above.

Computing a Range of Consecutive Terms

The goal is, given the polynomials \(P(x),Q(x)\) above, to compute \(\left\lbrack x^{\left\lbrack L,R\right)}\right\rbrack\dfrac{P(x)}{Q(x)}\). In the following computation, we only need to consider coefficients that "affect" the answer, which is key to the Bostan–Mori algorithm.

We can assume \(\deg{P}<\deg{Q}\); otherwise, we can reduce the problem to this case via one polynomial division with remainder.

First, consider a simpler problem:

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

We need to compute \(\left\lbrack x^{\left\lbrack L-\deg{Q},R\right)}\right\rbrack\dfrac{1}{Q(x)Q(-x)}\), then perform one multiplication and extract the coefficients of \(x^L,\dots ,x^{R-1}\). Let \(V(x^2)=Q(x)Q(-x)\). Then we only need to compute

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

to recover \(\left\lbrack x^{\left\lbrack L-\deg{Q},R\right)}\right\rbrack\dfrac{1}{Q(x)Q(-x)}\). Subsequently, we only need to compute \(\left\lbrack x^{\left\lbrack L-\deg{P},R\right)}\right\rbrack\dfrac{1}{Q(x)}\) and then multiply by \(P(x)\) once to obtain \(\left\lbrack x^{\left\lbrack L,R\right)}\right\rbrack\dfrac{P(x)}{Q(x)}\).

The above algorithm already works, but the time complexity of each recursion depends on \(R-L\). We hope to eliminate the dependence on \(R-L\) during the recursive computation. More specifically, first consider computing \(\left\lbrack x^{\left\lbrack L,L+\deg Q+1\right)}\right\rbrack \dfrac{1}{Q(x)}\). Consider

\[ \left\lbrack x^{\left\lbrack L,L+\deg Q+1\right)}\right\rbrack \frac{1}{Q(x)}=\left\lbrack x^{\left\lbrack L,L+\deg Q+1\right)}\right\rbrack \dfrac{1}{Q(x)Q(-x)}\cdot Q(-x) \]

We need to compute

\[ \left\lbrack x^{\left\lbrack L-\deg Q,L+\deg Q+1\right)}\right\rbrack \dfrac{1}{Q(x)Q(-x)} \]

Then for \(V(x^2)=Q(x)Q(-x)\), we only need to compute

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

This is because

\[ \left\lbrack x^{k}\right\rbrack\dfrac{1}{Q(x)Q(-x)}= \begin{cases} \left\lbrack x^{k/2}\right\rbrack\dfrac{1}{V(x)},&\text{if }k\equiv 0\pmod{2}, \\ 0,&\text{otherwise}. \end{cases} \]

We know \(L+\deg Q\) and \(L-\deg Q\) have the same parity, so

\[ \left\lceil \frac{L+\deg Q+1}{2}\right\rceil -\left\lceil \frac{L-\deg Q}{2}\right\rceil = \begin{cases} \deg Q+1,&\text{if }L+\deg Q\equiv 0\pmod{2}, \\ \deg Q,&\text{otherwise}. \end{cases} \]

Thus, we can write the pseudocode:

\[ \begin{array}{ll} &\textbf{Algorithm }\operatorname{Slice-Coefficients}(Q,L)\text{:} \\ &\textbf{Input}\text{: }Q(x)\in\mathbb{C}\left\lbrack x\right\rbrack,L\in\mathbb{Z}\text{.} \\ &\textbf{Output}\text{: }\left\lbrack x^{\left\lbrack L,L+\deg Q+1\right)}\right\rbrack Q(x)^{-1}\text{.} \\ 1&\textbf{if }L\leq 1\textbf{ then return }\left\lbrack x^{\left\lbrack L,L+\deg Q+1\right)}\right\rbrack Q(x)^{-1} \\ &\text{Use other algorithm to compute }Q(x)^{-1} \\ 2&V(x^2)\gets Q(x)Q(-x) \\ 3&k\gets \left\lceil \frac{L-\deg Q}{2}\right\rceil \\ 4&(t_k,\dots ,t_{k+\deg Q})\gets \operatorname{Slice-Coefficients}\left(V,k\right) \\ 5&T(x)\gets x^{(L-\deg Q)\bmod{2}}\sum_{j=0}^{\deg Q}t_{j+k}x^{2j} \\ 6&\textbf{return }\left\lbrack x^{\left\lbrack \deg Q,2\deg Q+1\right)}\right\rbrack T(x)Q(-x) \end{array} \]

However, this algorithm alone is insufficient; we need to find a new rational function representation and compute more coefficients.

Finding a New Rational Function Representation

We know \(Q(x)\) itself and a consecutive range of coefficients of \(Q(x)^{-1}\), e.g., \(\left\lbrack x^{\left\lbrack L,L+\deg Q\right)}\right\rbrack Q(x)^{-1}\) for \(L\geq 0\). We wish to compute \(\left\lbrack x^{\left\lbrack L+\deg Q,L+2\deg Q\right)}\right\rbrack Q(x)^{-1}\). This is equivalent to finding some \(P(x)\) with \(\deg P< \deg Q\) such that the first \(\deg Q\) terms of \(\dfrac{P(x)}{Q(x)}\) match \(\left\lbrack x^{\left\lbrack L,L+\deg Q\right)}\right\rbrack Q(x)^{-1}\). Simply put: the recurrence relation (the denominator of the rational function) remains unchanged; we are only changing the initial values (the numerator of the rational function).

Specifically, consider

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

Now we want to advance the recurrence by \(n\) terms, i.e.,

\[ \sum_{j\geq n}a_jx^{j-n}=\frac{P(x)}{Q(x)x^n}-\frac{Q(x)\sum_{j=0}^{n-1}a_jx^j}{Q(x)x^n} \]

We first use \(\operatorname{Slice-Coefficients}(Q,L-\deg{P})\) once to compute \(\left\lbrack x^{\left\lbrack L-\deg{P},L-\deg{P}+\deg{Q}+1\right)}\right\rbrack Q(x)^{-1}\). Then we extend and merge to obtain \(\left\lbrack x^{\left\lbrack L-\deg{P},L+\deg{Q}\right)}\right\rbrack Q(x)^{-1}\), and recompute a numerator such that

\[ \frac{\widetilde{P}(x)}{Q(x)}=\sum_{j\geq 0}\left(\left\lbrack x^{L+j}\right\rbrack \frac{P(x)}{Q(x)}\right)x^j \]

Finally, we use formal power series division to compute \(\left\lbrack x^{\left\lbrack 0,R-L\right)}\right\rbrack\dfrac{\widetilde{P}(x)}{Q(x)}\) in time \(O(\mathsf{M}(d)\log L+\mathsf{M}(R-L))\).

References

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

本页面最近更新：，更新历史
发现错误？想一起完善？在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用