# Beating binary powering for polynomial matrices

Alin Bostan  Inria and U. Paris-SaclayPalaiseauFrance<EMAIL_ADDRESS>,
Vincent Neiger  Sorbonne Université, CNRS, LIP6F-75005 Paris75252France
<EMAIL_ADDRESS>and  Sergey Yurkevich  University of Vienna and Inria
SaclayViennaAustria<EMAIL_ADDRESS>
(2023)

###### Abstract.

The $N$th power of a polynomial matrix of fixed size and degree can be
computed by binary powering as fast as multiplying two polynomials of linear
degree in $N$. When Fast Fourier Transform (FFT) is available, the resulting
arithmetic complexity is _softly linear_ in $N$, i.e. linear in $N$ with extra
logarithmic factors. We show that it is possible to beat binary powering, by
an algorithm whose complexity is _purely linear_ in $N$, even in absence of
FFT. The key result making this improvement possible is that the entries of
the $N$th power of a polynomial matrix satisfy linear differential equations
with polynomial coefficients whose orders and degrees are independent of $N$.
Similar algorithms are proposed for two related problems: computing the $N$th
term of a C-recursive sequence of polynomials, and modular exponentiation to
the power $N$ for bivariate polynomials.

Algebraic Algorithms; Computational Complexity; FFT; Binary Powering;
C-recursive Sequence; Rational Power Series; Linear Differential Equations;
Creative Telescoping; Polynomial Matrices.

The authors are supported by the ANR-19-CE40-0018 De Rerum Natura project and
by the joint ANR-FWF ANR-22-CE91-0007 EAGLES project. The third author was
supported by the ÖAW DOC fellowship P-26101.

††journalyear: 2023††copyright: acmlicensed††conference: Submitted; February
8; 2023††price: 15.00††ccs: Computing methodologies Algebraic algorithms

##  1\. Introduction

A sequence $(u_{n})_{n\geq 0}$ is called _C-recursive_ if it satisfies a
linear recurrence relation whose coefficients are constant with respect to
$n$. The famous sequence $(f_{n})_{n\geq 0}$ of Fibonacci numbers, defined by
the recurrence $f_{n+2}=f_{n+1}+f_{n}$ and the initial values
$f_{0}=0,f_{1}=1$, is perhaps the most basic example of a C-recursive sequence
after the geometric ones $(q^{n})_{n\geq 0}$. It is classical that the term
$f_{N}$ can be computed in $O(\log N)$ arithmetic operations, thus as fast as
$q^{N}$. This can be achieved by _binary powering_ for $q^{N}$, and in fact
for $f_{N}$ as well, since it is the top-right entry of $C^{N}$ where
$C=(\begin{smallmatrix}1&1\\\ 1&0\end{smallmatrix})$ is a $2\times 2$ matrix.
This idea generalizes to any C-recursive sequence $(u_{n})_{n\geq 0}$: a
recurrence of order $r\geq 1$ for $(u_{n})_{n\geq 0}$ can be encoded, via its
companion matrix, into an $r\times r$ matrix recurrence of order $1$. Then the
term $u_{N}$ of the sequence appears as an entry of the $N$th power of this
$r\times r$ companion matrix multiplied by the vector of initial values
$(u_{0},\dots,u_{r-1})$ (Miller and Brown, 1966; Fiduccia, 1985). Then $u_{N}$
can be computed in $O(\log N)$ arithmetic operations, and in $O(N\log(N))$ bit
operations if $(u_{n})_{n\geq 0}$ is an integer sequence, using fast integer
multiplication (Harvey and van der Hoeven, 2021). Here $r$ is considered as a
constant parameter, i.e.,
$r\in\mathchoice{O\left(1\right)}{O(1)}{O(1)}{O(1)}$.

Fibonacci polynomials $F_{n}(x)$ are a natural generalization of Fibonacci
numbers (see e.g. (Byrd, 1963)). They are defined by the recurrence

(1) $F_{n+2}(x)=xF_{n+1}(x)+F_{n}(x)\quad\text{for }n\geq 0$

and the initial values $F_{0}(x)=0,F_{1}(x)=1$. The first few terms are

$(F_{n})_{n\geq 0}=(0,1,x,x^{2}+1,x^{3}+2x,x^{4}+3x^{2}+1,\dots).$

Obviously, for all $n\geq 1$, the polynomial $F_{n}(x)$ is monic of degree
$n-1$ and the sum of its coefficients is $F_{n}(1)=f_{n}$.

Given $N\in\mathbb{N}$, the direct iterative algorithm for computing
$F_{N}(x)$ has arithmetic complexity $O(N^{2})$. It computes, for each $n\leq
N$, all the $n$ coefficients of the intermediate polynomial $F_{n}(x)$; in
total this amounts to $\Theta(N^{2})$ coefficients. Therefore, if one wants to
compute all of $(F_{0},\ldots,F_{N})$ then this direct method is optimal with
respect to the total arithmetic size of the output. However, it becomes
quadratic if one is only interested in determining $F_{N}(x)$ alone.

To compute the polynomial $F_{N}(x)$ faster, one can use, as in the scalar
case, the reformulation of the second-order recurrence (1) as a first-order
(polynomial) matrix recurrence:

(2) $\begin{pmatrix}F_{n+2}&F_{n+1}\\\
F_{n+1}&F_{n}\end{pmatrix}=\begin{pmatrix}x&1\\\
1&0\end{pmatrix}\begin{pmatrix}F_{n+1}&F_{n}\\\ F_{n}&F_{n-1}\end{pmatrix}.$

This shows that $F_{n}(x)$ is the top-right entry of the matrix $C(x)^{n}$,
where $C(x)$ is the $2\times 2$ polynomial matrix
$C(x)=(\begin{smallmatrix}x&1\\\ 1&0\end{smallmatrix})$. One can again compute
$C(x)^{N}$ using binary powering, whose costliest step is the multiplication
of two polynomial matrices each with degree about $N/2$. This leads to an
algorithm that finds $F_{N}(x)$ in arithmetic complexity
$O(\mathchoice{\mathsf{M}\left(N\right)}{\mathsf{M}(N)}{\mathsf{M}(N)}{\mathsf{M}(N)})$,
where
$\mathchoice{\mathsf{M}\left(N\right)}{\mathsf{M}(N)}{\mathsf{M}(N)}{\mathsf{M}(N)}$
denotes the arithmetic cost of polynomial multiplication in degree at most
$N$.

Using FFT-based polynomial multiplication, this amounts to a number of
arithmetic operations in the base field $\mathbb{K}$ which is quasi-linear in
$N$ (Cantor and Kaltofen, 1991). Not only does this compare favorably to the
complexity $O(N^{2})$ of the direct iterative algorithm, but this is even
quasi-optimal (i.e., optimal up to logarithmic factors) with respect to the
arithmetic size $\Theta(N)$ of the output polynomial $F_{N}(x)$.

In this context, the idea also generalizes to any C-recursive sequence
$(u_{n}(x))_{n\geq 0}$ of polynomials in $\mathbb{K}[x]$, which we will call
_polynomial C-recursive sequences_. Indeed, one can encode any recurrence of
arbitrary (but independent of $n$) order $r\geq 1$ and coefficients in
$\mathbb{K}[x]$ into a polynomial $r\times r$ matrix recurrence of order 1,
and the $N$th term of the sequence, $u_{N}(x)$, can be computed as an element
in the $N$th power of an $r\times r$ polynomial matrix multiplied by the
polynomial vector of initial values. Conversely, computing the $N$th power of
any polynomial matrix can be reduced to computing terms in polynomial
C-recursive sequences (see the introduction of Section 4). Binary powering
allows to solve both problems in
$O(\mathchoice{\mathsf{M}\left(N\right)}{\mathsf{M}(N)}{\mathsf{M}(N)}{\mathsf{M}(N)})$
arithmetic operations, and in $O(N^{2}\log(N))$ bit operations if
$\mathbb{K}=\mathbb{Q}$, considering both the recurrence order (or the matrix
size) $r$ and the recurrence degree (or the matrix degree) $d$ as constant
parameters, i.e., $r,d\in\mathchoice{O\left(1\right)}{O(1)}{O(1)}{O(1)}$. The
main question addressed in this article is:

> Can one achieve a better complexity for these tasks?

As far as scalar C-recursive sequences are concerned, the arithmetic
complexity $O(\log N)$ seems very difficult (if not impossible) to beat, but
it is perhaps not impossible to improve the bit complexity $O(N\log(N))$
towards $O(N)$. While we do not achieve this, our results provide polynomial
analogues for this type of improvement. As frequently noticed in computer
algebra, polynomials are “computationally easier” to deal with than integers.
In our case, philosophically, this comes from the fact that we can benefit
from an additional operation on polynomials: differentiation.

This possibly cryptic remark will hopefully become clear throughout the next
section. There, using the specific case of Fibonacci polynomials as a test
bench, we argue why it is indeed legitimate to hope for algorithms of
arithmetic complexity $O(N)$ for computing the $N$th term of a polynomial
C-recursive sequence.

#### Main result

Recall that a _C-finite sequence_ is a sequence $(u_{n})_{n\geq 0}$ of
elements $u_{n}$ in some ring $R$ which satisfies a recurrence equation

(3) $u_{n+r}=c_{r-1}u_{n+r-1}+\cdots+c_{0}u_{n}\quad\text{for all}\;n\geq 0,$

for $c_{0},\dots,c_{r-1}\in R$. In this work we consider _polynomial_ C-finite
sequences, i.e., the case $R=\mathbb{K}[x]$ for some (effective) field
$\mathbb{K}$ of characteristic zero; thus $u_{n}=u_{n}(x)\in\mathbb{K}[x]$.
The customary data structure for representing such a sequence consists of the
polynomials $c_{0}(x),\dots,c_{r-1}(x)$ defining the recurrence and the $r$
initial conditions $u_{0}(x),\dots,u_{r-1}(x)\in\mathbb{K}[x]$. The _order_ of
the recurrence is $r$ while its _degree_ is the maximum of the degrees of the
$c_{i}$’s.

###### Theorem 1.1.

Let $\mathbb{K}$ be an effective field of characteristic 0. Let $d$ and $r$ be
fixed positive integers. For each of the following problems, there exists an
algorithm solving it in $O(N)$ operations $(\pm,\times,\div)$ in $\mathbb{K}$:

  * SeqTerm: Given a polynomial C-finite sequence $(u_{n}(x))_{n\geq 0}$ of order and degree at most $r$ and $d$, compute the $N$th term $u_{N}(x)$.

  * BivModPow: Given polynomials $Q(x,y)$ and $P(x,y)$ in $\mathbb{K}[x,y]$ of degrees in $y$ and $x$ at most $r$ and $d$, with $P(x,y)$ monic in $y$, compute $Q(x,y)^{N}\bmod P(x,y)$.

  * PolMatPow: Given a square polynomial matrix $M(x)$ over $\mathbb{K}[x]$ of size and degree at most $r$ and $d$, compute $M(x)^{N}$.

Our algorithms for SeqTerm, BivModPow, and PolMatPow make essential use of
divisions in $\mathbb{K}$. We do not know if the complexity $O(N)$ can be
achieved using only operations $(+,-,\times)$ in $\mathbb{K}$.

#### Previous work

As already mentioned, the classical way of computing the $N$th term of a given
C-finite sequence goes via binary powering of the companion matrix, see e.g.
(Miller and Brown, 1966). Another algorithm, due to Fiduccia (Fiduccia, 1985),
works via binary powering in a polynomial quotient ring, and has a better
complexity with respect the order $r$, but not with respect to $N$. The
fastest existing algorithm (Bostan and Mori, 2021) for computing the $N$th
term of a C-finite sequence (counting arithmetic operations in $\mathbb{K}$)
beats Fiduccia’s algorithm by a constant factor. In our polynomial C-finite
case, and under the assumption
$r,d\in\mathchoice{O\left(1\right)}{O(1)}{O(1)}{O(1)}$, all these algorithms
achieve an arithmetic complexity in
$O(\mathchoice{\mathsf{M}\left(N\right)}{\mathsf{M}(N)}{\mathsf{M}(N)}{\mathsf{M}(N)})$.

Beyond this classical approach, the previous work on the aforementioned
problems consists of two distinct directions. The special case of Chebyshev
polynomials of the second kind $U_{n}(x)=(-i)^{n}F_{n+1}(2ix)$ (with
$F_{n}(x)$ the $n$th Fibonacci polynomial and $i$ the imaginary unit) was
considered in (Koepf, 1999) (and later in (Czirbusz, 2012)). These references
present various methods for the computation of the Chebyshev polynomials (of
the first and second kind) with arithmetic complexity ranging from $O(N)$ to
$O(N^{3})$. The results in (Koepf, 1999; Czirbusz, 2012) exploit the
particular structure of these polynomials; except for possibly other families
of classical orthogonal polynomials, for which explicit (hypergeometric)
formulas exist, the methods in (Koepf, 1999; Czirbusz, 2012) do not admit
obvious generalizations.

An idea closely connected to a fundamental building block of our algorithms is
explained in (Flajolet and Salvy, 1997, Pbm. 4). There, Flajolet and Salvy
exploit the fact that, given a polynomial $P(x)$ in $\mathbb{K}[x]$, the
coefficient sequence of the $n$th power $P(x)^{n}$ satisfies a linear
recurrence of order independent of $n$, and with coefficients in
$\mathbb{K}[x,n]$ of degree independent of $n$; this recurrence allows them to
compute (a selected coefficient of) $P(x)^{N}$ more efficiently than by binary
powering. This idea has been applied in (Bostan et al., 2007b, §8) to count
points on hyperelliptic curves over finite fields, with applications to
cryptography. The technique also yields a general solution to SeqTerm when
$r=1$.

#### Outline

The following observation generalizes that in (Flajolet and Salvy, 1997): the
coefficient sequence of the $n$th power of any _algebraic function_ satisfies
a recurrence of order and degree independent of $n$. From this, in Section 3,
we give algorithms for SeqTerm with cost
$\mathchoice{O\left(N\right)}{O(N)}{O(N)}{O(N)}$.

To complete the proof of Theorem 1.1, we design reductions between the three
problems. Obviously $\textsc{PolMatPow}\Rightarrow\textsc{SeqTerm}$, i.e., any
algorithm for PolMatPow with cost $O(N)$ induces one for SeqTerm with cost
$O(N)$ as well. Indeed, the $N$th term of a polynomial C-recursive sequence is
equal to an entry of the product of the $N$th power of a companion matrix and
the vector of initial values, and this polynomial matrix-vector multiplication
costs $O(N)$. Conversely, it also holds that
$\textsc{SeqTerm}\Rightarrow\textsc{PolMatPow}$. One natural way to see this
is to consider $r^{2}$ sequences corresponding to each entry of $M(x)^{n}$,
with recurrence given by the characteristic polynomial of $M(x)$; see the
introduction of Section 4. In Section 4.2, we give a more efficient algorithm
for this reduction, based on an algorithm for
$\textsc{SeqTerm}\Rightarrow\textsc{BivModPow}$ described in Section 4.1.

#### Basics of complexity and holonomic functions

Throughout the text, $\mathbb{K}$ denotes an effective field of characteristic
zero. For analyzing the performance of algorithms we will use the arithmetic
complexity model, meaning that arithmetic operations $(+,-,\times,\div)$ in
the base field $\mathbb{K}$ are counted at unit cost. As before,
$\mathchoice{\mathsf{M}\left(N\right)}{\mathsf{M}(N)}{\mathsf{M}(N)}{\mathsf{M}(N)}$
stands for the complexity of multiplication of two polynomials in
$\mathbb{K}[x]$ of degree at most $N$. With FFT-based multiplication
$\mathchoice{\mathsf{M}\left(N\right)}{\mathsf{M}(N)}{\mathsf{M}(N)}{\mathsf{M}(N)}\in
O(N\log(N)\log\log(N))$ (Cantor and Kaltofen, 1991), improved to $O(N\log(N))$
if $\mathbb{K}$ contains suitable roots of unity (Cooley and Tukey, 1965) or
if $\mathbb{K}$ is any finite field (Harvey and van der Hoeven, 2022).

A power series $f(x)\in\mathbb{K}[[x]]$ is called _D-finite_ (or _holonomic_)
if it satisfies a linear differential equation (ODE) of the form

(4) $q_{\ell}(x)f^{(\ell)}(x)+\cdots+q_{0}(x)f(x)=0,$

for some $q_{0}(x),\dots,q_{\ell}(x)\in\mathbb{K}[x]$ with $q_{\ell}(x)\neq
0$. Equivalently, writing $f(x)=\sum_{k\geq 0}f_{k}x^{k}$, the sequence
$(f_{k})_{k\geq 0}$ is _P-recursive_ , i.e., it satisfies a linear recurrence
equation (LRE)

$p_{s}(k)f_{k+s}+\cdots+p_{0}(k)f_{k}=0\quad\text{for all}\;k\geq 0,$

with polynomial coefficients $p_{0}(x),\dots,p_{s}(x)\in\mathbb{K}[x]$, and
$p_{s}\neq 0$. Note that $s$ and $\ell$ may differ in general, but
$s\leq\ell+\max_{i}(\deg q_{i}(x))$. It also holds that $\max_{i}(\deg
p_{i}(x))\leq\ell$.

It is often useful to write (4) as $Lf(x)=0$, where

$L=q_{\ell}(x)\partial_{x}^{\ell}+\cdots+q_{0}(x)$

is an element in the non-commutative Weyl algebra
$\mathbb{K}[x]\langle\partial_{x}\rangle$ of linear differential operators
with multiplication governed by the Leibniz rule
$\partial_{x}x=x\partial_{x}+1$. The _order_ $\ell$ of a differential operator
$L$ is the highest power of $\partial_{x}$ occurring in $L$, and the _degree_
of $L$ is the highest power of $x$. We recall that the least common left
multiple (LCLM) of two operators
$L_{1},L_{2}\in\mathbb{K}[x]\langle\partial_{x}\rangle$ is the unique monic
operator $L=\mathrm{LCLM}(L_{1},L_{2})$ of minimal order such that there exist
two nonzero operators $A,B\in\mathbb{K}(x)\langle\partial_{x}\rangle$ with
$L=AL_{1}=BL_{2}$. The LCLM can be computed efficiently (Bostan et al., 2012).

##  2\. The case of Fibonacci polynomials

Before solving the first part (SeqTerm) of Theorem 1.1 in general, we propose
in this section three different approaches that can be used to compute the
$N$th Fibonacci polynomial $F_{N}(x)$ in arithmetic complexity $O(N)$. Two of
these methods have the advantage that they generalize to the case of arbitrary
C-finite sequences.

The starting point of all that follows is the observation that the generating
function $F(x,y)\coloneqq\sum_{n\geq 0}F_{n}(x)y^{n}\in\mathbb{K}[x][[y]]$ of
the sequence $(F_{n}(x))_{n\geq 0}$ is _rational_ , and equal to
$y/(1-xy-y^{2})$.

###  2.1. First method via a closed-form expression

By using the partial fraction decomposition

$\frac{y}{1-xy-y^{2}}=\frac{1}{\varphi_{+}(x)-\varphi_{-}(x)}\cdot\left(\frac{1}{1-\varphi_{+}(x)y}-\frac{1}{1-\varphi_{-}(x)y}\right)$

where $\varphi_{\pm}(x)={(x\pm{\sqrt{x^{2}+4}})}/{2}$ are the roots of
$\varphi^{2}-x\varphi-1=0$, and by applying the geometric series, we get the
closed form expression

(5)
$F_{n}(x)=\frac{\varphi_{+}(x)^{n}-\varphi_{-}(x)^{n}}{\varphi_{+}(x)-\varphi_{-}(x)}\quad\text{for
all }n\geq 0.$

Now, using the binomial formula twice, we obtain the formula

(6) $F_{n}(x)=\frac{1}{2^{n-1}}\cdot\sum_{\ell\geq 0}4^{\ell}\left(\sum_{k\geq
0}\binom{n}{2k+1}\binom{k}{\ell}\right)x^{n-2\ell-1}.$

The identity (Gould, 1972, 3.121) implies a “magic” simplification:

(7) $\sum_{k\geq
0}\binom{n}{2k+1}\binom{k}{\ell}=2^{n-1-2\ell}\binom{n-\ell-1}{\ell}.$

In conclusion, from (6) and (7) it follows that

(8) $F_{n}(x)=\sum_{\ell\geq 0}\binom{n-\ell-1}{\ell}x^{n-2\ell-1}.$

With this expression at hand, it becomes transparent that one can compute
$F_{N}(x)$ efficiently. Indeed, by writing
$F_{N}(x)=\sum_{k=0}^{N-1}f_{k}x^{k}$, it follows from (8) that
$(f_{n})_{n\geq 0}$ satisfies the recurrence relation

(9) $f_{k+2}=\frac{(N+k+1)(N-k-1)}{4(k+1)(k+2)}f_{k}\quad\text{ for all }k\geq
0.$

Moreover, (8) also gives $(f_{0},f_{1})=(1,0)$ for odd $N$ and otherwise
$(f_{0},f_{1})=(0,N/2)$. With these initial conditions, it is now clear that
$F_{N}(x)$ can be computed in $O(N)$ by unrolling the recurrence (9).

As mentioned in the introduction, the analogue of formula (8) for the case of
Chebyshev polynomials of the first kind $T_{n}(x)$ was already exploited in
(Koepf, 1999, §1.9). The disadvantage of this approach is that for general
polynomial C-finite sequences there is no hope for a closed form expression
like (8).

###  2.2. Second method via algebraic substitution

There is another method for computing $F_{N}(x)$ in $O(N)$, which has the
advantage that it generalizes to any C-recursive sequence, as we will show in
Section 3.1.

The crucial remark (Lemma 3.2) is that since $\varphi_{\pm}(x)$ is algebraic,
$\varphi_{\pm}(x)^{n}$ satisfies a “small” linear differential equation: an
ODE of order and degree independent of $n$. The same holds true for
$1/(\varphi_{+}(x)-\varphi_{-}(x))$, therefore for $F_{n}(x)$ as well. More
precisely, $\varphi_{\pm}(x)^{n}$ satisfies the linear differential equation

$(x^{2}+4)y^{\prime\prime}(x)+xy^{\prime}(x)-n^{2}y(x)=0,$

and $1/(\varphi_{+}(x)-\varphi_{-}(x))=(x^{2}+4)^{-1/2}$ satisfies the ODE

$(x^{2}+4)y^{\prime}(x)+xy(x)=0.$

Using (5), it then follows that the polynomial $F_{n}(x)$ satisfies

(10) $(x^{2}+4)y^{\prime\prime}(x)+3xy^{\prime}(x)+(1-n^{2})y(x)=0.$

Writing $F_{N}(x)=\sum_{k=0}^{N-1}f_{k}x^{k}$, plugging into (10) for $n=N$
and extracting the $(k+2)$nd coefficient, it now follows that the sequence
$(f_{k})_{k\geq 0}$ satisfies recurrence (9). The initial conditions
$f_{0},f_{1}$ are given by $F_{N}(x)\bmod x^{2}$ which can be found in
complexity $O(\log N)$ by computing the $N$th power of the companion matrix
(2) in $\mathbb{K}[x]/(x^{2})$ by binary powering and reducing mod $x^{2}$ in
each step. As before, unrolling recurrence (9), with these initial terms,
provides a way to compute $F_{N}(x)$ in complexity $O(N)$.

###  2.3. Third method via Creative Telescoping

Writing $F(x,y)=y/(1-xy-y^{2})$ we are interested in a differential equation
for the coefficient of $y^{N}$ in $F(x,y)$. By Cauchy’s integral formula, we
have for sufficiently small $\epsilon>0$:

$F_{N}(x)=[y^{N}]F(x,y)=\frac{1}{2\pi
i}\oint_{|y|=\epsilon}\frac{y}{(1-xy-y^{2})y^{N+1}}\mathrm{d}y.$

Then the method of creative telescoping can be used to find a differential
equation for the integral above. For example, the command

        DEtools[Zeilberger](1/(1-x*y-y^2)/y^n, x, y, Dx);

in Maple immediately finds that

$\displaystyle\left((x^{2}+4)\partial_{x}^{2}+3x\partial_{x}+1-n^{2}\right)\frac{F(x,y)}{y^{n+1}}=\partial_{y}\left(\frac{F(x,y)}{y^{n}}C(x,y)\right),$

where $C(x,y)=\frac{n+1-nxy-(n-1)y^{2}}{1-xy-y^{2}}$. By Cauchy’s integral
theorem, the contour integral of the right-hand side vanishes, and thus (10)
follows. Then one can conclude in the same way as in the previous method and
compute $F_{N}(x)$ in complexity $O(N)$.

###  2.4. Comments on the three approaches

It is natural to ask ourselves what in these approaches was just luck, what
was truly specific to the particular example of the Fibonacci polynomials, and
what can be extended to the general case.

It is clear that the key for computing $F_{N}(x)$ in complexity $O(N)$ is the
existence of the recurrence (9) (or equivalently the ODE (10)). Even though
there is no hope for a closed form solution in general, we shall prove that
such a recurrence always exists for polynomial C-finite sequences. We should,
however, definitely be careful and avoid proving tautologic statements. Since
$u_{N}(x)$ is a polynomial, it does satisfy _a priori_ a first-order linear
differential equation with polynomial coefficients,
$u_{N}(x)y^{\prime}(x)-u_{N}^{\prime}(x)y(x)=0$, but this one is trivial for
our purposes. Indeed, converting this differential equation into a recurrence
satisfied by the sequence of coefficients of $u_{N}(x)$ yields a recurrence of
order $\deg(u_{N})$, which is obviously useless for computing the coefficients
of $u_{N}$. Therefore, we need to ensure existence of a recurrence/ODE whose
order and degree are independent of $N$. This is the purpose of the next
section. Specifically, in §3.1 we will explain how it can be found by
algebraic substitution (generalizing §2.2) and in §3.2 we will show that it
can also be found in general via creative telescoping (generalizing §2.3).

##  3\. Polynomial C-finite sequences

Recall that a polynomial C-finite sequence $(u_{n}(x))_{n\geq 0}$ is a
sequence of polynomials $u_{n}(x)\in\mathbb{K}[x]$ that satisfies a recurrence

(11) $u_{n+r}(x)=c_{r-1}(x)u_{n+r-1}(x)+\cdots+c_{0}(x)u_{n}(x),$

of some order $r\in\mathbb{N}$ and polynomial coefficients
$c_{0}(x),\dots,c_{r-1}(x)\in\mathbb{K}[x]$. The _degree_ of (11) is
$d=\max_{i}(\deg c_{i}(x))$. Eq. 11 defines the sequence $(u_{n}(x))_{n\geq
0}$ uniquely if $r$ initial terms $u_{0}(x),\dots,u_{r-1}(x)$ are prescribed.
The characteristic polynomial for (11) is defined as

$\chi(y)=y^{r}-c_{r-1}(x)y^{r-1}-\cdots-c_{1}(x)y-c_{0}(x)\in\mathbb{K}[x,y].$

It is well-known and not difficult to see that the generating function
$U(x,y)\coloneqq\sum_{n\geq 0}u_{n}(x)y^{n}$ is a rational function and given
by

(12) $U(x,y)=\frac{v_{0}(x)+\cdots+v_{r-1}(x)y^{r-1}}{y^{r}\chi(1/y)},$

with $v_{k}(x)\coloneqq u_{k}(x)-c_{r-1}(x)u_{k-1}(x)-\cdots-
c_{r-k}(x)u_{0}(x)$.

Let $a_{1}(x),\dots,a_{k}(x)\in\overline{\mathbb{K}(x)}$ be the roots of
$\chi(y)=0$, and $m_{1},\dots,m_{k}$ be their multiplicities. It follows from
the partial fraction decomposition and geometric series that any sequence
$(u_{n}(x))_{n\geq 0}$ satisfying (11) is of the form

(13) $u_{n}(x)=q_{1}(n,x)a_{1}(x)^{n}+\cdots+q_{k}(n,x)a_{k}(x)^{n},$

where $k\leq r$ and each $q_{i}(n,x)\in\mathbb{K}(a_{1}(x),\dots,a_{n}(x))[n]$
is a polynomial in $n$ of degree at most $m_{i}-1$, for $i=1,\dots,k$.

###  3.1. Computing $u_{N}(x)$ in $O(N)$

By generalizing the ideas of Section 2.2, it is not difficult to prove that
the $n$th term of a polynomial C-recursive sequence $(u_{n}(x))_{n\geq 0}$
satisfies a linear ODE whose degree and order are independent of $n$, and
consequently, that there exists a linear recurrence relation for the
coefficient sequence of $u_{n}(x)$ whose order (say $s$) and degree are again
independent of $n$. Then, for a given $N\in\mathbb{N}$, first computing
initial terms by binary powering of the companion matrix in
$\mathbb{K}[x]/(x^{s})$ and then unrolling this recurrence for $n=N$, we
achieve a complexity $O(N)$ for the computation of $u_{N}(x)$.

###### Theorem 3.1.

Let $(u_{n}(x))_{n\geq 0}$ be a polynomial C-recursive sequence. Then there
exists $L_{n}\in\mathbb{K}[n,x]\langle\partial_{x}\rangle$ with order and
degree independent of $n$, and such that $L_{n}(u_{n}(x))=0$. Consequently,
writing $u_{n}(x)=\sum_{k\geq 0}c_{n,k}x^{k}$, there exist, for some
$s\in\mathbb{N}$ independent of $n$, polynomials
$p_{0}(n,x),\dots,p_{s}(n,x)\in\mathbb{K}[n,x]$ of degrees independent of $n$,
and such that the sequence $(c_{n,k})_{k\geq 0}$ satisfies the recurrence

(14) $p_{s}(n,k)c_{n,k+s}+\cdots+p_{0}(n,k)c_{n,k}=0,\quad k\geq 0.$

In the theorem above it is crucial that neither the order nor the degree of
$L_{n}$ depend on $n$. Since each $u_{n}(x)$ is a polynomial, it is a
tautology to say that it satisfies _some_ linear differential equation with
polynomial coefficients: one may simply take $L=\partial_{x}^{\alpha}$, where
$\alpha>\deg(u_{n}(x))$ or $L=u_{n}(x)\partial_{x}-u_{n}^{\prime}(x)$.
However, it is a nontrivial fact that $u_{n}(x)$ satisfies an ODE of the form

$p_{\ell}(n,x)u_{n}^{(\ell)}(x)+\dots+p_{0}(n,x)u_{n}(x)=0$

for some $p_{i}(n,x)\in\mathbb{K}[n,x]$ (with $\ell$ and $\deg_{x}p_{i}$
independent of $n$).

The most direct proof of Theorem 3.1 uses the explicit expression (13) for
$u_{n}(x)$ and the following classical fact about algebraic substitution into
D-finite functions. Recall that a function $a(x)$ is called _algebraic_ over
$\mathbb{K}(x)$ if it satisfies a non-trivial polynomial relation
$P(x,a(x))=0$ for some $P(x,y)\in\mathbb{K}[x,y]$. Size and complexity bounds
on differential equations for algebraic functions, and more generally on
algebraic substitution, are given in (Bostan et al., 2007a; Kauers and
Pogudin, 2017).

Algorithm 1 SeqTermAS$((u_{n})_{n},N)$

1:A polynomial C-finite sequence $(u_{n}(x))_{n\geq 0}$ given by (11) with
initial conditions, and $N\in\mathbb{N}$.

2:The polynomial $u_{N}(x)$.

3:$\chi(y)\leftarrow$ the characteristic polynomial of $(u_{n}(x))_{n\geq 0}$

4:$a_{1}(x),\ldots,a_{k}(x)\leftarrow$ the roots of $\chi(y)$

5:Compute minimal polynomials for $q_{1}(x,n),\dots,q_{k}(x,n)$
$\in\mathbb{K}(a_{1}(x),\dots,a_{k}(x))[n]$ such that (13) holds.

6:For each $i$ deduce an ODE
$L_{i,n}\in\mathbb{K}[n,x]\langle\partial_{x}\rangle$ with order and degree
independent of $n$ such that $L_{i,n}(q_{i}(x,n)a_{i}(x)^{n})=0$

7:$L_{n}\leftarrow\mathrm{LCLM}(L_{1,n},\dots,L_{k,n})\in\mathbb{K}[n,x]\langle\partial_{x}\rangle$

8:Compute a recurrence $p_{s}(n,k)c_{n,k+s}+\dots+p_{0}(n,k)c_{n,k}=0$
satisfied by any solution $f_{n}(x)=\sum_{k\geq 0}c_{n,k}x^{k}$ of $L_{n}y=0$

9:Using binary powering of the companion matrix of the initial recurrence mod
$x^{s}$, compute the values $c_{N,0},\dots,c_{N,s-1}$

10:Unroll the recurrence from 8 for $n=N$ and with initial terms from 9

11:return $\sum_{k=0}^{Nd}c_{N,k}x^{k}$, where $d=\deg_{x}(\chi(y))$

###### Lemma 3.2.

Let $a(x)$ be an algebraic function over $\mathbb{K}(x)$ and let $g(x)$ be
D-finite. Then $f(x)=g(a(x))$ is D-finite. In particular, $a(x)^{n}$ satisfies
a linear ODE of order and degree independent of $n$.

###### Proof.

The first part is a classical result, see for example (Stanley, 1980, Thm.
2.7). In the proof one shows that the vector space spanned over
$\mathbb{K}(x)$ by $(f^{(i)}(x))_{i\geq 0}$ is finite-dimensional over
$\mathbb{K}(x,a(x))$ which is itself finite-dimensional over $\mathbb{K}(x)$.
For the second part, it is enough to set $g(x)=x^{n}$ which satisfies
$xg^{\prime}(x)=ng(x)$. ∎

###### Example 3.3.

Like in Section 2 let $\varphi_{\pm}(x)=(x\pm\sqrt{x^{2}+4})/2$ be the roots
of $y(x)^{2}+xy(x)-1=0$. Then $\varphi_{\pm}(x)^{n}$ satisfy the ODE

$(x^{2}+4)y^{\prime\prime}(x)+xy^{\prime}(x)-n^{2}y(x)=0.$

###### Proof of Theorem 3.1.

Write $u_{n}(x)$ as in (13). By Lemma 3.2, each $a_{i}(x)^{n}$ satisfies a
linear differential equation of order and degree independent of $n$, hence the
same holds for $q_{i}(n,x)a_{i}(x)^{n}$, and finally for $u_{n}(x)$. It
follows that the coefficient sequence of $u_{n}(x)$ is P-recursive with order
and degree independent of $n$. ∎

Since all steps in the proofs above are effective and independent of $N$, this
leads to Algorithm 1. Its 3 to 8 can be seen as “precomputations” since they
do not depend on $N$. As already mentioned, 9 has complexity $O(\log N)$ and
10 has complexity $O(N)$. Thus, Algorithm 1 solves SeqTerm in complexity
$O(N)$, up to a potential issue during the unrolling at 10 of the recurrence
from 8. Indeed, this unrolling may be impossible for some values $k$, namely
those for which $p_{s}(N,k)$ vanishes. We will explain how to overcome this
problem in Section 3.3.

For practical applications, however, computing the polynomials $q_{i}(x,n)$ in
5 as well as the LCLM in 7 is algorithmically somewhat cumbersome. Thus,
generalizing the approach in Section 2.3, we now propose a variant of
Algorithm 1 which replaces 3 to 7 by an algorithm based on creative
telescoping.

###  3.2. Computing $L_{n}$ with Creative Telescoping

Let $U(x,y)=\sum_{n\geq 0}u_{n}(x)y^{n}\in\mathbb{K}[x][[y]]$ be the
generating function (12) of $(u_{n}(x))_{n\geq 0}$. The sequence is C-finite,
so $U(x,y)$ is a rational function. Moreover, the Cauchy integral formula
implies that

(15) $u_{n}(x)=\frac{1}{2\pi
i}\oint_{|y|=\epsilon}\frac{U(x,y)}{y^{n+1}}\mathrm{d}y.$

A _telescoper_ of $U(x,y)/y^{n+1}$ is a differential operator

$L=p_{k}(x)\partial_{x}^{k}+\cdots+p_{0}(x)\in\mathbb{K}[x]\langle\partial_{x}\rangle,$

such that $L$ applied to $U(x,y)/y^{n+1}$ is $\partial_{y}(C(x,y))$ for some
rational function $C(x,y)$ called the _certificate_. By the Cauchy integral
theorem, $\oint_{|y|=\epsilon}\partial_{y}(C(x,y))\mathrm{d}y=0$, and it
follows that $Lu_{n}(x)=0$, i.e., $L$ yields a differential equation for
$u_{n}(x)$. In this section we will prove that for $U(x,y)/y^{n+1}$ there
exists a telescoper $L_{n}\in\mathbb{K}[n,x]\langle\partial_{x}\rangle$ whose
order and degree do not depend on $n$. Our proof relies on reduction-based
creative telescoping and repeatedly uses the _Hermite reduction_ algorithm
(Bostan et al., 2010, 2013a, 2018).

We now introduce the necessary definitions and recall the Hermite reduction
method. For a more detailed introduction, a full complexity analysis, and
applications of reduction-based creative telescoping to integration of
bivariate rational functions, we refer to (Bostan et al., 2010). Let
$\mathbb{L}=\mathbb{K}(x)$. For a polynomial $Q(y)\in\mathbb{L}[y]$, let
$Q=Q_{1}Q_{2}^{2}\cdots Q_{k}^{k}$ be its squarefree factorization and let
$Q^{*}=Q_{1}\cdots Q_{k}$ denote the squarefree part of $Q$. We set
$Q^{-}\coloneqq Q/Q^{*}$. Recall that, given $P,Q\in\mathbb{L}[y]$, the
_Hermite reduction_ algorithm computes two polynomials $A,a\in\mathbb{L}[y]$
with $\deg_{y}a<\deg_{y}Q^{*}$ such that

$\frac{P}{Q}=\partial_{y}\left(\frac{A}{Q^{-}}\right)+\frac{a}{Q^{*}}.$

Given a bivariate rational function $H(x,y)=P(y)/Q(y)\in\mathbb{K}(x,y)$, one
may compute the Hermite reduction $(A_{i},a_{i})$ of $\partial_{x}^{i}H$ for
$i=0,1,\dots$. Since $\deg_{y}a_{i}$ is uniformly bounded by
$d^{*}=\deg_{y}Q^{*}$ for each $i$, the $d^{*}+1$ functions
$\\{a_{i}(x,y)\colon 0\leq i\leq d^{*}\\}$ will be linearly dependent over
$\mathbb{K}(x)$. Hence one can find
$q_{0}(x),\cdots,q_{d^{*}}(x)\in\mathbb{K}(x)$ not all zero, such that
$\sum_{i=0}^{d^{*}}q_{i}(x)a_{i}(x)=0$. It follows then that
$L=\sum_{i=0}^{d^{*}}q_{i}(x)\partial_{x}^{i}$ is a telescoper for $H$.

This procedure cannot be directly applied to $U(x,y)/y^{n+1}$ if $n$ is an
indeterminate. At the same time, if $n=N\in\mathbb{N}$ is fixed, it is a
priori not obvious that $\deg_{x}q_{i}(x)$ will be independent of $N$.
Moreover, the complexity of the algorithm will depend on $N$, which we want to
avoid. As we will now explain, to achieve this, one should see
$U(x,y)/y^{n+1}$ not as a rational function in $x$ and $y$ with potentially
large degree in the numerator, but as a hyperexponential function with the
parameter $n$ appearing solely as a coefficient in the logarithmic derivative.
Recall that $H(x,y)$ is called _hyperexponential_ if $\partial_{x}H/H$ and
$\partial_{y}H/H$ belong to $\mathbb{K}(x,y)$.

For hyperexponential functions, the Almkvist-Zeilberger algorithm (Almkvist
and Zeilberger, 1990) was the first practical method to find telescopers and
certificates. Indeed, as we mentioned in Section 2.3, the command

        DEtools[Zeilberger](1/(1-x*y-y^2)/y^n, x, y, Dx);

in Maple immediately finds the differential equation for the $n$th Fibonacci
polynomial for a variable $n$. Note that if $n$ is specialized to an integer
$N$ before the execution of the command above, the implemented algorithm
becomes slower as $N$ growths.

It is, however, not clear that the Almkvist-Zeilberger algorithm applied to
$U(x,y)/y^{n+1}$ will always find a telescoper whose degree and order are
independent of $n$, even though we know from Section 3.1 that an ODE with this
property exists. Therefore, to have a complete algorithm based on creative
telescoping, we will invoke the reduction-based method for hyperexponential
functions first introduced and analyzed in (Bostan et al., 2013a). Using the
implementation of the latter work, the command in Maple

        HermiteTelescoping(1/(1-x*y-y^2)/y^n, x, y, Dx);

also immediately finds the correct ODE for $F_{n}(x)$. The practical advantage
for our purpose of using the reduction based algorithm in comparison to the
Almkvist-Zeilberger method is shown in Section 5 (Table 1). The theoretical
advantage comes from the following lemma, which guarantees that the algorithm
will find a telescoper for $U(x,y)/y^{n+1}$, and consequently an ODE for
$u_{n}(x)$, whose order and degree do not depend on $n$.

###### Lemma 3.4.

Let $P(y)\in\mathbb{L}[n,y]$ and $Q(y)\in\mathbb{L}[y]$ with $Q(0)\neq 0$. Set
$d_{n}\coloneqq\deg_{n}P(y)$, $d^{*}\coloneqq\deg_{y}Q^{*}(y)$ and let $k$ be
the highest pure power in the square free factorization of $Q(y)$. Then there
exist $B(n,y)\in\mathbb{L}(n)[y]$ and $b(n,y)\in\mathbb{L}[n,y]$ with
$\deg_{y}b(n,y)\leq d^{*}$ and $\deg_{n}b(n,y)\leq d_{n}+k$ such that

(16)
$\frac{P(y)}{Q(y)y^{n+1}}=\partial_{y}\left(\frac{B(n,y)}{Q^{-}(y)y^{n}}\right)+\frac{b(n,y)}{Q^{*}(y)y^{n+1}}.$

###### Proof.

We are going prove the statement by induction on
$d^{-}\coloneqq\deg_{y}Q^{-}$. If $d^{-}=0$, then $Q^{*}=Q$ and the Euclidean
division gives $P=P_{1}Q+b_{1}$ with $\deg_{y}b_{1}<d^{*}$. Moreover,

$\frac{P_{1}(y)}{y^{n+1}}=\partial_{y}\left(\frac{B_{1}(n,y)}{y^{n}}\right),$

where $B_{1}(n,y)$ is $P_{1}(y)$ with the $k$th coefficient $p_{k}$ replaced
by $p_{k}/(k-n)$. Setting $B=B_{1}$ and $b=b_{1}$ proves the induction basis.

Now assume that $d^{-}>0$ and note that

$\partial_{y}\left(\frac{B(n,y)}{Q^{-}(y)y^{n}}\right)=y^{-n}\partial_{y}\left(\frac{B(n,y)}{Q^{-}(y)}\right)-y^{-n-1}n\frac{B(n,y)}{Q^{-}(y)},$

so equation (16) is equivalent to

(17)
$\frac{P(y)}{Q(y)y}=\partial_{y}\left(\frac{B(n,y)}{Q^{-}(y)}\right)-n\frac{B(n,y)}{Q^{-}(y)y}+\frac{b(y)}{Q^{*}(y)y}.$

The Hermite reduction algorithm applied to $\frac{P(y)}{Q(y)y}$ yields
polynomials $A(y),a(y)$ with $\deg_{y}a(y)\leq d^{*}$ such that

(18)
$\frac{B(y)}{Q(y)y}=\partial_{y}\left(\frac{A(y)}{Q^{-}(y)}\right)+\frac{a(y)}{Q^{*}(y)y}.$

Comparing (17) and (18), we now look at

$H(y)\coloneqq\frac{a(y)}{Q^{*}(y)y}+n\frac{A(y)}{Q^{-}(y)y}.$

The denominator of $H(y)$ is $R(y)y\coloneqq\mathrm{lcm}(Q^{*},Q^{-})y$.
Clearly, $R^{*}=Q^{*}$ and $\deg_{y}R^{-}<d^{-}$. The highest pure power in
the square free factorization of $R(y)$ is at most $k-1$ and the degree of the
numerator in $n$ of $H(y)$ is at most $d_{n}+1$. Hence, by induction, we may
write

$\frac{a(y)}{Q^{*}(y)y}+n\frac{A(y)}{Q^{-}(y)y}=\partial_{y}\left(\frac{C(n,y)}{Q^{-}(y)}\right)-n\frac{C(n,y)}{Q^{-}(y)y}+\frac{c(n,y)}{Q^{*}(y)y}$

with $\deg_{y}c(n,y)<d^{*}$ and $\deg_{n}c(n,y)\leq d_{n}+k$. Setting
$B(n,y)=A(y)+C(n,y)$ and $b(y)=a(y)+c(n,y)$ finishes the proof. ∎

The proof of Lemma 3.4 induces an algorithm for the computation of $B(n,y)$
and $b(n,y)$ given $P(y),Q(y)\in\mathbb{K}[x,y]$ such that (16) holds,
$\deg_{y}b\leq\deg_{y}Q^{*}$ and also $\deg_{n}b$ bounded in terms of $Q$. It
can be seen as a special case of the procedure HermiteReduction in (Bostan et
al., 2013a). The ODE for $u_{n}(x)$ can be now found as in Algorithm 2.

Algorithm 2 TelescNthTerm$(U(x,y))$

1:A rational function $U(x,y)\in\mathbb{K}(x,y)\cap\mathbb{K}[[x,y]]$.

2:A diff. operator in $\mathbb{K}[n,x]\langle\partial_{x}\rangle$ for
$u_{n}(x)=[y^{n}]U(x,y)$.

3:Write $U(x,y)=P(x,y)/Q(x,y)$ and let $d^{*}=\deg_{y}Q^{*}(x,y)$

4:For each $i=0,\dots,d^{*}$ compute the polynomial $b_{i}(n,x,y)=b(n,y)$ as
in Lemma 3.4 applied to $\partial_{x}^{i}U(x,y)/y^{n+1}$

5:Find a linear relation of $\\{b_{i}(n,x,y)\colon 0\leq i\leq d^{*}\\}$ over
$\mathbb{K}(n,x)$, that is polynomials $q_{0}(n,x),\dots,q_{d^{*}}(n,x)$ not
all zero with $\sum_{i=0}^{d^{*}}q_{i}(n,x)b_{i}(n,x,y)=0$

6:Return the differential operator
$\sum_{i=0}^{d^{*}}q_{i}(n,x)\partial_{x}^{i}\in\mathbb{K}[n,x]\langle\partial_{x}\rangle$

Note that, as in the usual reduction based creative telecoping, the linear
relation at 5 exists because $\deg_{y}(b_{i}(n,x,y))$ is uniformly bounded by
$d^{*}$. Writing $U(x,y)=P(x,y)/Q(x,y)$, the operator
$L=\sum_{i=0}^{d^{*}}q_{i}(n,x)\partial_{x}^{i}$ annihilates
$u_{n}(x)=[y^{n}]U(x,y)$ since

$\displaystyle 2\pi i\cdot L_{n}u_{n}(x)=L_{n}\oint\frac{U(x,y)}{y^{n+1}}$
$\displaystyle\mathrm{d}y=\oint L_{n}\frac{P(x,y)}{Q(x,y)y^{n+1}}\mathrm{d}y$
$\displaystyle=\oint\partial_{y}\frac{\sum_{i=1}^{d^{*}}q_{i}(x,y)B_{i}(n,y)}{Q^{-}(x,y)y^{n}}$
$\displaystyle\mathrm{d}y+\oint\frac{\sum_{i=0}^{d^{*}}q_{i}(n,x)b_{i}(n,x,y)}{Q^{*}(x,y)y^{n+1}}\mathrm{d}y;$

the first integral vanishes by Cauchy’s integral theorem, and the second
integral vanishes by construction of the $q_{i}(n,x)$.

This provides a variant for 3 to 7 of Algorithm 1, as described in Algorithm
3.

Algorithm 3 SeqTermCT$((u_{n})_{n},N)$

1:A polynomial C-finite sequence $(u_{n}(x))_{n\geq 0}$ given by (11) with
initial conditions, and $N\in\mathbb{N}$.

2:The polynomial $u_{N}(x)$.

3:$U(x,y)\leftarrow$ the rational generating function of $u_{n}(x)$ in (12)

4:$L_{n}\leftarrow$ $\textsc{TelescNthTerm}(U(x,y))$

5:$\triangleright$ _follow 8 to 11 of Algorithm 1_

As above, a potential issue is that the unrolling step is only possible for
$p_{s}(N,k)\neq 0$. The next section deals with this problem.

###  3.3. The singular case

In this section we discuss the potential issue of our algorithm that can occur
if the sequence for the coefficients of $u_{n}(x)$ cannot be unrolled due to
singularities. We shall first highlight this problem and its solution by means
of an example.

Consider the polynomial C-recurrent sequence $u_{n}(x)$ given by

$u_{n+3}(x)-(x^{2}+x+2)u_{n+2}(x)+x(x^{2}+2x+2)u_{n+1}(x)-2x^{3}u_{n}(x)=0,$

for all $n\geq 0$ with initial conditions
$u_{0}=3,u_{1}=x^{2}+x+2,u_{2}=x^{4}+x^{2}+4$. The characteristic polynomial
of the defining recurrence is easily computed and turns out to factor
completely over $\mathbb{K}[x][y]$:

$\chi(x,y)=(y-2)(y-x)(y-x^{2}).$

With the initial conditions and after a partial fraction decomposition it
follows that the generating function of $u_{n}(x)$ is given by

$U(x,y)=\frac{1}{1-2y}+\frac{1}{1-x^{2}y}+\frac{1}{1-xy}.$

Hence, the solution is $u_{n}(x)=2^{n}+x^{n}+x^{2n}$, and can be even written
down in $O(\log(N))$ arithmetic operations for any $N$. However, as we shall
explain now, the direct application of any of the methods described earlier
fails.

According to Theorem 3.1, $u_{n}(x)$ satisfies an ODE whose degree and order
are independent of $n$. Indeed, using creating telescoping one quickly finds
an annihilator for $u_{n}(x)=\oint U(x,y)/y^{n+1}\mathrm{d}x$:

$\bigl{(}x^{2}\partial_{x}^{3}-3x(n-1)\partial_{x}^{2}+(2n-1)(n-1)\partial_{x}\bigr{)}u_{n}(x)=0.$

Converting this ODE to a recurrence for the coefficient sequence of
$u_{n}(x)=\sum_{k\geq 0}c_{n,k}x^{k}$ we find

(19) $(2n-k)(n-k)kc_{n,k}=0,\quad k\geq 0.$

In other words, $c_{n,k}=0$ for all $k\in\mathbb{N}$ except
$k\in\\{0,n,2n\\}$. In order to “unroll” this recurrence we need to know
$c_{n,0},c_{n,n}$ and $c_{n,2n}$. However, it is not immediately clear how to
compute those terms for $n=N$ in $O(N)$ arithmetic operations from the initial
input (without using the explicit solution).

We propose the following easily generalizable solution: consider
$v_{n}(x)=u_{n}(x+1)$. Then the ODE for $v_{n}(x)$ is given by

$\bigl{(}(x+1)^{2}\partial_{x}^{3}-3(x+1)(n-1)\partial_{x}^{2}+(2n-1)(n-1)\partial_{x}\bigr{)}v_{n}(x)=0,$

and for the coefficient sequence of $v_{n}(x)=\sum_{k\geq 0}d_{n,k}x^{k}$ we
find

$(k+1)(k+2)d_{n,k+2}-(k+1)(3N-2k-1)d_{n,k+1}+(2n-k)(n-k)d_{n,k}=0.$

Now the leading coefficient of the recurrence is $(k+1)(k+2)\neq 0$, so we can
easily unroll it after determining the first two terms, by computing them via
binary powering of the corresponding companion matrix mod $x^{2}$. Having
computed $v_{N}(x)$, it remains to find $u_{N}(x)=v_{N}(x-1)$. Note that
expanding the polynomial results in an
$O(\mathchoice{\mathsf{M}\left(N\right)}{\mathsf{M}(N)}{\mathsf{M}(N)}{\mathsf{M}(N)})$
algorithm. However, recall from (19) that we only need to compute $c_{N,N}$
and $c_{N,2N}$, or, in other words, the coefficients of $x^{N}$ and $x^{2N}$
in $v_{N}(x-1)$. For any $i$ it holds that

(20) $c_{N,i}=\sum_{k\geq 0}d_{N,k}\binom{k}{i}(-1)^{k-i},$

and the sum is finite because $v_{N}(x)$ is a polynomial. Clearly, it can be
computed in complexity $O(N)$ for any $i$.

Generally speaking, an issue with unrolling the recurrence for
$(c_{n,k})_{k\geq 0}$ occurs if the roots of the leading polynomial are
positive integers that depend on $n$. Let $S$ be the set of these roots. Note
that the size of $S$ is independent of $n$ and $S$ can be non-empty only if
the ODE for $u_{n}(x)$ is singular at 0 (that is, if $q_{\ell}(x)$ in (4)
vanishes at $0$). In this case, one can always define $v_{n}(x)=u_{n}(x+c)$
for $c\in\mathbb{K}$ a non-singular point of the ODE ($q_{\ell}(c)\neq 0$).
Then the coefficients $d_{n,k}$ of $v_{n}(x)$ can be computed from $O(1)$
initial conditions via unrolling a recurrence. Using the formula (20) (with
$-c$ instead $-1$) and the fact that $v_{n}(x)$ is a polynomial, one can
compute the coefficients $c_{N,i}$ for $i\in S$. With these, it is then
possible to unroll the recurrence for $(c_{N,k})_{k\geq 0}$ and find
$u_{N}(x)$ in complexity $O(N)$.

##  4\. Impact on polynomial matrix power

Here is an algorithm for PolMatPow using SeqTerm. Let $M(x)$ in
$\mathbb{K}[x]_{\leq d}^{r\times r}$ and $p_{i,j,n}(x)$ be the $(i,j)$ entry
of $M(x)^{n}$, for $n\geq 0$ and $i$ and $j$ in $\\{1,\ldots,r\\}$. The
sequence $(p_{i,j,n}(x))_{n\geq 0}$ is polynomial C-recursive, with a
recurrence given by the characteristic polynomial

$\chi_{M}(x,y)\coloneqq\det(y\,I_{r}-M(x))=y^{r}-c_{r-1}(x)y^{r-1}-\cdots-
c_{0}(x).$

That is,
$p_{i,j,n+r}(x)=c_{r-1}(x)p_{i,j,n+r-1}(x)+\cdots+c_{0}(x)p_{i,j,n}(x)$ for
all $n\geq 0$. Thus, to compute $M(x)^{N}$, it is enough to find
$\chi_{M}(x,y)$ (in $O(1)$, i.e. independent of $N$), to compute the
polynomials $p_{i,j,n}(x)$ for $1\leq i,j\leq r$ and $0\leq n<r$ (also in
$O(1)$) and to return the entries $p_{i,j,N}(x)$ of $M(x)^{N}$ via SeqTerm. As
such, this approach uses $r^{2}$ calls to SeqTerm, with total cost
$\mathchoice{O\left(N\right)}{O(N)}{O(N)}{O(N)}$.

This section describes an algorithm for PolMatPow which uses only $r$ such
calls, through a direct reduction to BivModPow (see Section 4.2). Our solution
for BivModPow, via $r$ calls to SeqTerm, is presented in Section 4.1 and
completes the proof of Theorem 1.1.

###  4.1. Computing bivariate modular powers

Let $\mathbb{L}=\mathbb{K}[x]$ and $P,Q\in\mathbb{L}[y]$. Assume that $P$,
seen as a univariate polynomial in $y$ of degree $r$, is monic. For
$N\in\mathbb{N}$, Euclidean division in $\mathbb{L}[y]$ ensures the existence
of unique $S,R\in\mathbb{L}[y]$ such that $\deg_{y}(R)<r$ and $Q^{N}=SP+R$.
The polynomial $R$ is $Q^{N}\bmod P$. Assume that $P$ and $Q$ are fixed, and
let $d\coloneqq\deg_{x}(P)$ (which is thus in $O(1)$). Then, writing
$R=\sum_{i=0}^{d-1}r_{i}(x)y^{i}$, it holds that $\deg_{x}r_{i}(x)=O(N)$. The
efficient computation of $R$ when $Q=y$, given $P(x,y)$ and $N$, is the first
step for proving BivModPow in Theorem 1.1.

We shall first illustrate the connection of SeqTerm and BivModPow by means of
an example. Let $P(x,y)=y^{2}-xy-1$ and $Q(x,y)=y$, i.e., we are looking for
$F_{n-1}(x),F_{n}(x)\in\mathbb{L}$ such that

(21) $y^{n}=S(x,y)(y^{2}-xy-1)+yF_{n}(x)+F_{n-1}(x),$

for some polynomial $S(x,y)\in\mathbb{L}[y]$. Replace $y$ by $1/y$ in (21) and
then multiply by $y^{n+1}/(1-xy-y^{2})$ to obtain

$\frac{y}{1-xy-y^{2}}=Q(x,1/y)y^{n-1}+y^{n}\frac{F_{n}(x)+yF_{n-1}(x)}{1-xy-y^{2}}.$

Now observe that $\deg_{y}(Q(x,1/y)y^{n-1})\leq n-1$, hence by extracting the
$n$th and $(n+1)$st coefficients,

$\displaystyle F_{n}(x)$ $\displaystyle=[y^{n}]\frac{y}{1-xy-y^{2}}\text{
and}$ $\displaystyle F_{n-1}(x)+xF_{n}(x)$
$\displaystyle=[y^{n+1}]\frac{y}{1-xy-y^{2}}.$

We conclude that $F_{k}(x)$ is the $k$th Fibonacci polynomial, for $k=n$ and
$k=n-1$. In particular, each $F_{k}(x)$ satisfies a linear recurrence with
constant polynomials and can be found in $O(k)$ by SeqTerm.

This strategy, outlined on an example, generalizes in the obvious way.
Explicitly, we have the following lemma (see (Bostan and Mori, 2021, Lem. 2)).

###### Lemma 4.1.

Let $P\in\mathbb{K}[x,y]$ and $r\coloneqq\deg_{y}(P)$, with $P(x,0)\neq 0$ and
reversal $\bar{P}(x,y)\coloneqq y^{r}P(x,\frac{1}{y})$. Write
$\frac{1}{\bar{P}(x,y)}\eqqcolon\sum_{k\geq 0}u_{k}(x)y^{k}.$ Finally, let
$v(x,y)=(u_{N-r+1}(x)+\cdots+u_{N}(x)y^{r-1})\bar{P}(x,y)\bmod y^{r}$. Then
$y^{N}\bmod P(x,y)=v(1/y)y^{r-1}$.

The sequence $(u_{k}(x))_{k\geq 0}$ in Lemma 4.1 is C-recursive because its
generating function is rational. Hence, using SeqTerm, the $r=O(1)$ many terms
$u_{N-r+1}(x),\dots,u_{N}(x)$ can be computed in complexity $O(N)$. It follows
that the case $Q(x,y)=y$ of BivModPow can be solved in $O(N)$ steps as well.

Finally, the computation of $Q(x,y)^{N}\bmod P(x,y)$ can be reduced to
$y^{N}\bmod P(x,y)$ with a resultant precomputation (see Lemma 4.2). This
leads to Algorithm 4, which solves BivModPow in
$\mathchoice{O\left(N\right)}{O(N)}{O(N)}{O(N)}$.

###### Lemma 4.2.

Let $P(y),Q(y)\in\mathbb{L}[y]$. Define $A(t),B(t)\in\mathbb{L}[t]$ by
$A(t)=\mathrm{Res}_{y}(P(y),t-Q(y))$ and $B(t)=t^{N}\bmod A(t)$. Then

$Q(y)^{N}\bmod P(y)=B(Q(y))\bmod P(y).$

###### Proof.

By the definition of the resultant, $A(t)=\prod_{i}(t-Q(a_{i}))$ where
$a_{i}\in\overline{\mathbb{L}}$ are the solutions of $P(y)=0$. Hence, $P(y)$
divides $A(Q(y))$, which, by construction, divides $B(Q(y))-Q(y)^{N}$. ∎

Algorithm 4 BivModPow$(P(x,y),Q(x,y),N)$

1:$P(x,y),Q(x,y)\in\mathbb{L}[y]$ with $P(x,y)$ monic in $y$, and
$N\in\mathbb{N}$.

2:$Q(x,y)^{N}\bmod P(x,y)$.

3:$A(x,t)\leftarrow\mathrm{Res}_{y}(P(x,y),t-Q(x,y))$

4:$\bar{A}(x,t)\leftarrow t^{r}A(x,1/t)$, where $r\coloneqq\deg_{t}A(x,t)$

5:for $i=N-r+1,\dots,N$ do

6: $u_{i}(x)\leftarrow[t^{i}]\frac{1}{\bar{A}(x,t)}$ using SeqTerm

7:$u(x,t)\leftarrow u_{N-r+1}(x)+\cdots+u_{N}(x)t^{r-1}$

8:$v(x,t)\leftarrow u(x,t)\bar{A}(x,t)\bmod t^{r}$; $B(x,t)\leftarrow
v(1/t)t^{r-1}$

9:return $B(x,Q(x,y))\bmod P(x,y)$

###  4.2. Computing polynomial matrix powers

Let $M\in\mathbb{K}[x]^{r\times r}$ be an $r\times r$ polynomial matrix of
degree at most $d$. Its power $M^{N}$ has degree at most $Nd\in O(N)$. Let
$P(x,y)$ be the characteristic polynomial of $M$. Since $P(x,M)=0$ by the
Cayley-Hamilton theorem, we get $M^{N}=R(x,M)$ where $R(x,y)=y^{N}\bmod
P(x,y)$. The polynomial $R$ can be computed in $O(N)$ via BivModPow. Then
evaluating $R(x,y)$ at $y=M(x)$ has cost $O(N)$ because $\deg_{x}(R)\in O(N)$,
$\deg_{y}(R)=r\in O(1)$ and $\deg_{x}(M)=d\in O(1)$. Hence Algorithm 5 is
correct and has complexity $O(N)$.

Algorithm 5 PolMatPow$(M,N)$

1:matrix $M(x)\in\mathbb{K}[x]^{r\times r}$, integer $N\in\mathbb{N}$.

2:$M(x)^{N}\in\mathbb{K}[x]^{r\times r}$.

3:$P(x,y)\leftarrow$ the characteristic polynomial of $M(x)$

4:$R(x,y)\leftarrow y^{N}\bmod P(x,y)$ $\triangleright$ _instance of
BivModPow_

5:return $R(x,M(x))$

##  5\. Experiments

The main precomputation step for all our algorithms consists in starting with
a rational function $U\in\mathbb{K}(x,y)\cap\mathbb{K}[[x,y]]$ and in finding
a differential operator $L_{n}$ that annihilates $u_{n}(x)=[y^{n}]U(x,y)$ and
whose degree and order are independent of $n$. For this task, in practice, we
may either use the method described in §3.1, or creative telescoping
algorithms for hyperexponential functions. Table 1 summarizes timings for a
variety of implementations.

The table reveals that, among these implementations, the fastest one for
computing a telescoper of $U(x,y)/y^{n+1}$ is the reduction-based creative
telescoping in Maple. More specifically, redct is the fastest, followed by HT.
The implementation in _ore_algebra_ (Kauers and Mezzarobba, 2019) in SageMath
competes best with reduction-based methods.

Table 2 gives timings of an efficient implementation of the remaining stages
after precomputations: computing initial terms (IT), and unrolling (UR). We
observe that IT takes negligible time compared to UR, except for extreme
parameter ranges where, simultaneously, $r$ and $d$ are large and $N$ is
small; this is expected since these ranges correspond to cases where the order
of the recurrence to be unrolled is close to $N$. We also see that binary
powering is always slower, often by a factor more than $5$, than the addition
of IT and UR. The speed-up factor is summarized in Figs. 1, 2 and 3; as
expected it grows when $N$ grows, with $r$ and $d$ fixed.

For large $N$, in most of the reported cases, performing both the
precomputation and IT$+$UR is much faster than using binary powering. Still,
this is not always true, e.g. for $r=5$. One has to keep in mind that redct is
not implemented in low-level Maple, and targets rational coefficients: for a
more meaningful assessment of the precomputation part, it would be interesting
to have an implementation of creative telescoping which is fully optimized and
specialized to coefficients in a word-size prime field.

##  6\. Conclusion and future work

We have shown that it is possible to beat, both in theory and in practice, the
very basic and powerful binary powering method for computing: (i) powers of
polynomial matrices, (ii) terms in polynomial C-finite sequences and (iii)
modular exponentiation for bivariate polynomials. We describe below several
lines of work, including possible optimizations and generalizations, that we
leave for future investigations.

#### More detailed complexity analysis.

The first and most natural direction for future work is to analyze and improve
the complexity of the algorithms in Theorem 1.1 with respect to the parameters
$r$ and $d$. For simplicity, these parameters were assumed to be $O(1)$ in
this work. For the $N$th power of a polynomial matrix $M(x)$ of size $r$ and
degree $d$, binary splitting has arithmetic complexity
$O(\mathchoice{\mathsf{M}\left(Nd\right)}{\mathsf{M}(Nd)}{\mathsf{M}(Nd)}{\mathsf{M}(Nd)}r^{2}+{Nd}r^{\omega})$,
where $\omega\in[2,3]$ is a feasible exponent of matrix multiplication over
$\mathbb{K}$. With our approach, it is legitimate to target a differential
equation satisfied by the entries of $M(x)^{N}$ of order $r$ with coefficients
in $x$ of degree $O(dr^{3})$, yielding a recurrence of order $O(dr^{3})$ and
coefficients in $n$ of degree at most $r$. For large $N$, the complexity of
the new algorithm would thus be
$O(Ndr^{2}\mathchoice{\mathsf{M}\left(r\right)}{\mathsf{M}(r)}{\mathsf{M}(r)}{\mathsf{M}(r)})$.
Using different ODEs, of order $O(r)$ and coefficients of degree $O(dr^{2})$
could even lead to a complexity
$O(Ndr\mathchoice{\mathsf{M}\left(r\right)}{\mathsf{M}(r)}{\mathsf{M}(r)}{\mathsf{M}(r)})$.

#### The $K$th coefficient of the $N$th term

For some (large) integers $N,K\in\mathbb{N}$, one might be interested in
computing the single coefficient $[x^{K}y^{N}]U(x,y)$ of a rational function
$U\in\mathbb{K}(x,y)\cap\mathbb{K}[[x,y]]$. Equivalently it is natural to
wonder: How fast can one compute the $K$th coefficient of the $N$th term of a
C-recursive sequence $(u_{n}(x))_{n\geq 0}$? Using our method, a recurrence
with initial conditions for the coefficients of $u_{N}(x)$ can be deduced in
$O(\log N)$ arithmetic operations. Then (assuming that the recurrence is non-
singular) the $K$th coefficient can be found in
$O(\mathchoice{\mathsf{M}\left(\sqrt{K}\right)}{\mathsf{M}(\sqrt{K})}{\mathsf{M}(\sqrt{K})}{\mathsf{M}(\sqrt{K})})$
operations by using baby-steps/giant-steps (Chudnovsky and Chudnovsky, 1988;
Bostan et al., 2007b). We expect that, at least under a genericity assumption,
this problem can be solved in arithmetic complexity
$O(\log(N)+\mathchoice{\mathsf{M}\left(\sqrt{K}\right)}{\mathsf{M}(\sqrt{K})}{\mathsf{M}(\sqrt{K})}{\mathsf{M}(\sqrt{K})})$
which is a big improvement compared to the previous best $O(N+K)$ by (Massazza
and Radicioni, 2005).

#### Polynomial P-recursive sequences

A somewhat related task is to study the analogous problem to SeqTerm for
polynomial P-recursive sequences, that is for
$(u_{n}(x))_{n}\in\mathbb{K}[x]^{\mathbb{N}}$ satisfying

$p_{r}(x,n)u_{n+r}(x)+\cdots+p_{0}(x,n)u_{n}(x)=0,$

for $p_{i}(x,n)\in\mathbb{K}[x,n]$. We expect that, at least under a
genericity assumption, a generalization of Lemma 3.4 (based on results in
(Bostan et al., 2018)) should exist, implying in particular that $u_{N}(x)$
satisfies a linear ODE of order and degree independent of $N$. Generalizing
this even further, one might study the Creative Telescoping problem for
rational functions of the form
$H(\textbf{x})=\frac{P(x_{1},\dots,x_{s})}{Q(x_{1},\dots,x_{s})R(x_{1},\dots,x_{s})^{n}}.$
We expect that (at least generically) the minimal telecoper for
$H(\textbf{x})$ has order and degree independent of $n$ and can be found via a
Griffiths-Dwork reduction type approach, based on ideas from (Bostan et al.,
2013b).

#### Connection to the Jordan–Chevalley decomposition

A different approach for computing powers of matrices uses the
_Jordan–Chevalley decomposition_ (also called _SN decomposition_), see e.g.
(Hsieh et al., 1996; Elaydi and Harris, 1998; Schmidt, 2000; Couty et al.,
2011). It ensures that any polynomial matrix $M\in\mathbb{K}[x]^{r\times r}$
can be written as $M=S+Z,$ where $S\in\mathbb{K}(x)^{r\times r}$ is
diagonalizable over $\overline{\mathbb{K}(x)}$, $Z\in\mathbb{K}(x)^{r\times
r}$ is nilpotent, and $SZ=ZS$. From this decomposition it follows that
$M^{N}=\sum_{i=0}^{r-1}\binom{N}{i}S^{N-i}Z^{i}$. After a change of basis,
this reduces to computing a power of a diagonal matrix with algebraic
functions coefficients. Using Lemma 3.2 this can be performed efficiently in
$O(N)$ operations. It would be certainly interesting to compare this approach
with the other methods.

#### A PDE approach for SeqTerm

There is yet another method to deduce recurrence (9). The starting point is
that the generating function $F(x,y)=y/(1-xy-y^{2})$ of $F_{n}(x)$ satisfies
the linear PDE

(22) $(x^{2}+4)\frac{\partial^{2}F}{\partial x^{2}}+3x\frac{\partial
F}{\partial x}-y^{2}\frac{\partial^{2}F}{\partial y^{2}}-y\frac{\partial
F}{\partial y}+F=0,$

and extracting the coefficient of $x^{k}y^{n}$ in (22) immediately gives (9).
More generally, such a PDE translates into a recurrence if it is linear with
polynomial coefficients in $x$ and $y$ and if additionally all terms of the
form $x^{i}y^{\ell}\frac{\partial^{k}F}{\partial
x^{k}}\frac{\partial^{j}F}{\partial y^{j}}$ have $\ell=j$. A dimension
counting argument in the spirit of (Lipshitz, 1988, Lem. 3) proves that such a
PDE exists for _any_ rational function $F(x,y)$. The existence proof is
effective and amounts to linear algebra. A natural question is whether it is
possible to compute such a PDE via Creative Telescoping (either Almkvist-
Zeilberger (Almkvist and Zeilberger, 1990) or reduction based (Bostan et al.,
2010)), and how the corresponding method compares to the aforementioned ones.

Table 1. Timings in seconds for creative telescoping to find a telescoper
$L_{n}$ of $P(x,y)/(y^{n+1}Q(x,y))$. Here $P(x,y)/Q(x,y)$ is the generating
function for the sequence of the top-right entry of the powers of a randomly
chosen matrix in $\mathbb{F}_{p}[x]^{r\times r}$ of degree $d$, for a 50-bit
prime $p$, with $Q(x,y)$ the $y$-reversal of the characteristic polynomial of
this matrix. The order of $L_{n}$ is $\ell$, its degree in $n$ is
$\mathrm{d}_{n}$, and $\mathrm{d}_{x}=\deg_{x}(L_{n})$. A blank space means
that the computation took more than 1000 seconds. We observe empirically that
the degree in $x$ is $dr(r+1)(2r-1)/2-r(r-1)$ while its degree in $n$ is
$(r-1)(r+2)/2$; this is expected asymptotically by (Bostan et al., 2010, Thm.
25) and Lemma 3.4, because $\deg_{y}Q(x,y)=r$ and $\deg_{x}Q(x,y)=dr$. The
tested implementations are
• in Maple: redct (Bostan et al., 2018); HermiteTelescoping (HT) (Bostan et
al., 2013b); Zeilberger (ZB) (Almkvist and Zeilberger, 1990) in DEtools;
creative_telescoping (c_t) (Chyzak, 2000);
• in SageMath: creative telescoping (ct) from the _ore_algebra_ package
(Kauers and Mezzarobba, 2019);
• In Mathematica: FindCreativeTelescoping (FCT), CreativeTelescoping (CT), and HermiteTelescoping (HCT), see (Koutschan, 2010).  |  | Maple | Sage | Mathematica | $\ell$ | $\mathrm{d}_{n}$ | $\mathrm{d}_{x}$
---|---|---|---|---|---|---|---
$r$ | $d$ | redct | HT | ZB | c_t | ct | FCT | CT | HCT |  |  |
| 2 | 0.0 | 0.1 | 0.0 | 0.1 | 0.5 | 0.2 | 0.2 | 0.2 | 2 | 2 | 16
2 | 4 | 0.0 | 0.0 | 0.0 | 0.1 | 0.6 | 0.4 | 0.4 | 0.3 | 2 | 2 | 34
| 6 | 0.0 | 0.0 | 0.0 | 0.1 | 0.6 | 0.7 | 0.5 | 0.5 | 2 | 2 | 52
| 8 | 0.0 | 0.0 | 0.0 | 0.1 | 0.8 | 1.0 | 0.7 | 0.7 | 2 | 2 | 70
| 1 | 0.0 | 0.2 | 0.0 | 0.5 | 2.0 | 2.0 | 1.3 | 1.3 | 3 | 5 | 24
| 2 | 0.0 | 0.1 | 0.8 | 3.4 | 3.1 | 4.0 | 2.6 | 2.5 | 3 | 5 | 54
3 | 3 | 0.1 | 0.2 | 0.8 | 9.3 | 5.6 | 10 | 5.7 | 5.4 | 3 | 5 | 84
| 4 | 0.1 | 0.5 | 18 | 19 | 8.2 | 17 | 9.4 | 8.9 | 3 | 5 | 114
| 5 | 0.2 | 1.1 | 5.1 | 32 | 12 | 25 | 14 | 14 | 3 | 5 | 144
| 6 | 0.5 | 1.7 | 9.8 | 49 | 17 | 35 | 19 | 20 | 3 | 5 | 174
| 1 | 0.4 | 2.9 | 23 | 117 | 20 | 31 | 25 | 25 | 4 | 9 | 58
| 2 | 1.7 | 17 | 410 | 749 | 45 | 101 | 96 | 95 | 4 | 9 | 128
4 | 3 | 4.4 | 43 |  |  | 89 | 295 | 376 | 373 | 4 | 9 | 198
| 4 | 12 | 82 |  |  | 172 | 388 | 752 | 693 | 4 | 9 | 268
| 5 | 18 | 128 |  |  | 280 | 635 |  |  | 4 | 9 | 338
| 1 | 11 | 34 | 538 |  | 163 | 847 | 780 |  | 5 | 14 | 115
5 | 2 | 64 | 183 |  |  | 515 |  |  |  | 5 | 14 | 250
| 3 | 159 | 526 |  |  |  |  |  |  | 5 | 14 | 385
| 4 | 345 |  |  |  |  |  |  |  | 5 | 14 | 520
Table 2. Timings in seconds, using the C++ library NTL (Shoup, 2021) and PML
(Hyun et al., 2019), for computing the top-right entry of the $N$th power of a
randomly chosen matrix in $\mathbb{F}_{p}[x]^{r\times r}$ of degree $d$, for a
50-bit prime $p$. The first tested method is to directly apply binary powering
(BP); in the present context, the polynomial matrix multiplication of PML is
based on evaluation-interpolation and 3-prime FFT. The second tested method
uses Algorithm 3 and we do not count “precomputations” (already showed in
Table 1), i.e. we only report timings for the two non-negligible steps that
depend on $N$, namely 10 (UR, unrolling) and 9 (IT, initial terms) from
Algorithm 1.

|  | $N=2^{10}$ | $N=2^{12}$ | $N=2^{14}$ | $N=2^{16}$ | $N=2^{18}$ | $N=2^{20}$ | $N=2^{22}$
---|---|---|---|---|---|---|---|---
$r$ | $d$ | BP | UR | IT | BP | UR | IT | BP | UR | IT | BP | UR | IT | BP | UR | IT | BP | UR | IT | BP | UR | IT
| 2 | 1.2e-3 | 5.7e-4 | 3.7e-5 | 5.3e-3 | 2.4e-3 | 4.3e-5 | 2.5e-2 | 9.7e-3 | 4.9e-5 | 1.1e-1 | 3.9e-2 | 5.5e-5 | 5.3e-1 | 1.5e-1 | 6.2e-5 | 3.3e+0 | 6.2e-1 | 6.7e-5 | 1.5e+1 | 2.5e+0 | 7.5e-5
2 | 4 | 2.6e-3 | 1.3e-3 | 7.8e-5 | 1.2e-2 | 5.2e-3 | 9.4e-5 | 5.2e-2 | 2.1e-2 | 1.1e-4 | 2.4e-1 | 8.4e-2 | 1.3e-4 | 1.4e+0 | 3.4e-1 | 1.4e-4 | 7.2e+0 | 1.4e+0 | 1.6e-4 | 3.1e+1 | 5.4e+0 | 1.8e-4
| 6 | 3.8e-3 | 2.1e-3 | 1.2e-4 | 1.7e-2 | 8.7e-3 | 1.5e-4 | 7.9e-2 | 3.5e-2 | 1.8e-4 | 3.6e-1 | 1.4e-1 | 2.1e-4 | 2.3e+0 | 5.5e-1 | 2.4e-4 | 1.0e+1 | 2.2e+0 | 2.7e-4 | 4.6e+1 | 8.9e+0 | 3.0e-4
| 8 | 5.3e-3 | 3.1e-3 | 1.9e-4 | 2.4e-2 | 1.2e-2 | 2.4e-4 | 1.1e-1 | 5.0e-2 | 2.8e-4 | 5.3e-1 | 2.0e-1 | 3.3e-4 | 3.3e+0 | 8.0e-1 | 3.8e-4 | 1.5e+1 | 3.2e+0 | 4.3e-4 | 7.0e+1 | 1.2e+1 | 4.9e-4
| 1 | 1.4e-3 | 3.0e-4 | 1.3e-4 | 6.0e-3 | 1.3e-3 | 1.7e-4 | 2.6e-2 | 5.5e-3 | 2.1e-4 | 1.2e-1 | 2.2e-2 | 2.4e-4 | 5.8e-1 | 8.8e-2 | 2.8e-4 | 3.4e+0 | 3.5e-1 | 3.1e-4 | 1.6e+1 | 1.4e+0 | 3.5e-4
| 2 | 2.9e-3 | 7.8e-4 | 4.0e-4 | 1.2e-2 | 3.2e-3 | 5.3e-4 | 5.6e-2 | 1.3e-2 | 6.5e-4 | 2.6e-1 | 5.2e-2 | 7.8e-4 | 1.5e+0 | 2.1e-1 | 9.1e-4 | 7.6e+0 | 8.4e-1 | 1.0e-3 | 3.4e+1 | 3.3e+0 | 1.2e-3
3 | 3 | 4.3e-3 | 1.4e-3 | 7.4e-4 | 1.9e-2 | 5.8e-3 | 9.9e-4 | 8.4e-2 | 2.3e-2 | 1.2e-3 | 3.9e-1 | 9.3e-2 | 1.5e-3 | 2.2e+0 | 3.7e-1 | 1.7e-3 | 1.1e+1 | 1.5e+0 | 2.0e-3 | 4.9e+1 | 6.0e+0 | 2.2e-3
| 4 | 6.0e-3 | 2.1e-3 | 8.0e-4 | 2.6e-2 | 8.8e-3 | 1.0e-3 | 1.2e-1 | 3.5e-2 | 1.3e-3 | 5.8e-1 | 1.4e-1 | 1.5e-3 | 3.5e+0 | 5.7e-1 | 1.8e-3 | 1.7e+1 | 2.3e+0 | 2.0e-3 | 7.1e+1 | 9.1e+0 | 2.3e-3
| 5 | 7.4e-3 | 3.0e-3 | 1.0e-3 | 3.3e-2 | 1.2e-2 | 1.3e-3 | 1.5e-1 | 5.0e-2 | 1.7e-3 | 7.2e-1 | 2.0e-1 | 2.0e-3 | 4.3e+0 | 7.9e-1 | 2.3e-3 | 2.0e+1 | 3.2e+0 | 2.6e-3 | 8.8e+1 | 1.3e+1 | 2.9e-3
| 6 | 9.1e-3 | 4.0e-3 | 1.2e-3 | 4.0e-2 | 1.6e-2 | 1.6e-3 | 1.8e-1 | 6.6e-2 | 1.9e-3 | 8.2e-1 | 2.7e-1 | 2.3e-3 | 5.3e+0 | 1.1e+0 | 2.7e-3 | 2.3e+1 | 4.2e+0 | 3.1e-3 | 1.1e+2 | 1.7e+1 | 3.4e-3
| 1 | 2.7e-3 | 4.2e-4 | 7.8e-4 | 1.1e-2 | 1.8e-3 | 1.1e-3 | 4.9e-2 | 7.5e-3 | 1.4e-3 | 2.2e-1 | 3.0e-2 | 1.7e-3 | 1.1e+0 | 1.2e-1 | 2.0e-3 | 6.2e+0 | 4.8e-1 | 2.3e-3 | 2.9e+1 | 1.9e+0 | 2.6e-3
| 2 | 5.5e-3 | 1.2e-3 | 1.3e-3 | 2.4e-2 | 5.2e-3 | 1.8e-3 | 1.1e-1 | 2.1e-2 | 2.3e-3 | 4.9e-1 | 8.6e-2 | 2.8e-3 | 2.8e+0 | 3.4e-1 | 3.2e-3 | 1.4e+1 | 1.4e+0 | 3.7e-3 | 6.2e+1 | 5.5e+0 | 4.2e-3
4 | 3 | 8.2e-3 | 2.4e-3 | 2.1e-3 | 3.6e-2 | 1.0e-2 | 2.9e-3 | 1.6e-1 | 4.2e-2 | 3.7e-3 | 7.3e-1 | 1.7e-1 | 4.5e-3 | 4.4e+0 | 6.7e-1 | 5.3e-3 | 2.1e+1 | 2.7e+0 | 6.1e-3 | 9.3e+1 | 1.1e+1 | 6.9e-3
| 4 | 1.1e-2 | 4.1e-3 | 3.0e-3 | 5.0e-2 | 1.7e-2 | 4.1e-3 | 2.3e-1 | 6.9e-2 | 5.2e-3 | 1.1e+0 | 2.8e-1 | 6.4e-3 | 6.6e+0 | 1.1e+0 | 7.5e-3 | 3.1e+1 | 4.5e+0 | 8.6e-3 | 1.3e+2 | 1.8e+1 | 9.7e-3
| 5 | 1.4e-2 | 6.0e-3 | 3.7e-3 | 6.3e-2 | 2.5e-2 | 5.2e-3 | 2.8e-1 | 1.0e-1 | 6.6e-3 | 1.3e+0 | 4.1e-1 | 8.0e-3 | 7.7e+0 | 1.6e+0 | 9.4e-3 | 3.6e+1 | 6.5e+0 | 1.1e-2 | 1.6e+2 | 2.6e+1 | 1.2e-2
| 1 | 4.4e-3 | 6.0e-4 | 1.8e-3 | 1.8e-2 | 2.7e-3 | 2.5e-3 | 8.2e-2 | 1.1e-2 | 3.3e-3 | 3.7e-1 | 4.5e-2 | 4.1e-3 | 1.7e+0 | 1.8e-1 | 4.9e-3 | 1.0e+1 | 7.3e-1 | 5.7e-3 | 4.7e+1 | 2.9e+0 | 6.5e-3
5 | 2 | 9.1e-3 | 2.0e-3 | 3.5e-3 | 3.9e-2 | 9.1e-3 | 5.1e-3 | 1.8e-1 | 3.7e-2 | 6.7e-3 | 8.1e-1 | 1.5e-1 | 8.3e-3 | 4.6e+0 | 6.0e-1 | 9.9e-3 | 2.3e+1 | 2.4e+0 | 1.2e-2 | 1.0e+2 | 9.6e+0 | 1.3e-2
| 3 | 1.3e-2 | 4.3e-3 | 5.7e-3 | 5.8e-2 | 1.9e-2 | 8.3e-3 | 2.6e-1 | 7.8e-2 | 1.1e-2 | 1.2e+0 | 3.2e-1 | 1.4e-2 | 7.1e+0 | 1.3e+0 | 1.6e-2 | 3.4e+1 | 5.1e+0 | 1.9e-2 | 1.5e+2 | 2.0e+1 | 2.2e-2
| 4 | 1.8e-2 | 7.4e-3 | 7.8e-3 | 8.0e-2 | 3.3e-2 | 1.2e-2 | 3.8e-1 | 1.3e-1 | 1.5e-2 | 1.8e+0 | 5.4e-1 | 1.9e-2 | 1.1e+1 | 2.2e+0 | 2.3e-2 | 4.9e+1 | 8.7e+0 | 2.6e-2 | 2.1e+2 | 3.5e+1 | 3.0e-2

Figure 1. Speed-up versus binary powering, not counting precomputations, for
$r=2\ldots 7$, $N=2^{8},2^{10},\ldots,2^{22}$, and fixed $d=2$. Figure 2.
Speed-up versus binary powering, not counting precomputations, for $d=2\ldots
7$, $N=2^{8},2^{10},\ldots,2^{22}$, and fixed $r=4$. Figure 3. Speed-up versus
binary powering, not counting precomputations, for $r=2\ldots 7$, $d=1\ldots
7$, and fixed $N=2^{22}$.

## References

  * (1)
  * Almkvist and Zeilberger (1990) Gert Almkvist and Doron Zeilberger. 1990.  The method of differentiating under the integral sign.  _J. Symbolic Comput._ 10, 6 (1990), 571–591.  https://doi.org/10.1016/S0747-7171(08)80159-9
  * Bostan et al. (2010) Alin Bostan, Shaoshi Chen, Frédéric Chyzak, and Ziming Li. 2010.  Complexity of creative telescoping for bivariate rational functions. In _ISSAC’10_. ACM, 203–210.  https://doi.org/10.1145/1837934.1837975
  * Bostan et al. (2013a) Alin Bostan, Shaoshi Chen, Frédéric Chyzak, Ziming Li, and Guoce Xin. 2013a.  Hermite reduction and creative telescoping for hyperexponential functions. In _ISSAC’13_. ACM, 77–84.  https://doi.org/10.1145/2465506.2465946
  * Bostan et al. (2018) Alin Bostan, Frédéric Chyzak, Pierre Lairez, and Bruno Salvy. 2018.  Generalized Hermite reduction, creative telescoping and definite integration of D-finite functions. In _ISSAC’18_. ACM, 95–102.  https://doi.org/10.1145/3208976.3208992
  * Bostan et al. (2007a) Alin Bostan, Frédéric Chyzak, Grégoire Lecerf, Bruno Salvy, and Éric Schost. 2007a.  Differential equations for algebraic functions.  In _ISSAC’07_. ACM, 25–32.  https://doi.org/10.1145/1277548.1277553
  * Bostan et al. (2012) Alin Bostan, Frédéric Chyzak, Ziming Li, and Bruno Salvy. 2012.  Fast computation of common left multiples of linear ordinary differential operators. In _ISSAC’12_. ACM, 99–106.  https://doi.org/10.1145/2442829.2442847
  * Bostan et al. (2007b) Alin Bostan, Pierrick Gaudry, and Éric Schost. 2007b.  Linear recurrences with polynomial coefficients and application to integer factorization and Cartier-Manin operator.  _SIAM J. Comput._ 36, 6 (2007), 1777–1806.  https://doi.org/10.1137/S0097539704443793
  * Bostan et al. (2013b) Alin Bostan, Pierre Lairez, and Bruno Salvy. 2013b.  Creative telescoping for rational functions using the Griffiths-Dwork method. In _ISSAC’13_. ACM, 93–100.  https://doi.org/10.1145/2465506.2465935
  * Bostan and Mori (2021) Alin Bostan and Ryuhei Mori. 2021.  A Simple and Fast Algorithm for Computing the $N$-th Term of a Linearly Recurrent Sequence. In _SOSA’21 (Symposium on Simplicity in Algorithms)_. SIAM.  https://doi.org/10.1137/1.9781611976496.14
  * Byrd (1963) Paul F. Byrd. 1963\.  Expansion of analytic functions in polynomials associated with Fibonacci numbers.  _Fibonacci Quart._ 1, 1 (1963), 16–29.
  * Cantor and Kaltofen (1991) David G. Cantor and Erich Kaltofen. 1991.  On fast multiplication of polynomials over arbitrary algebras.  _Acta Inform._ 28, 7 (1991), 693–701.  https://doi.org/10.1007/BF01178683
  * Chudnovsky and Chudnovsky (1988) D. V. Chudnovsky and G. V. Chudnovsky. 1988.  Approximations and complex multiplication according to Ramanujan.  In _Ramanujan revisited (Urbana-Champaign, Ill., 1987)_. Academic Press, Boston, MA, 375–472.
  * Chyzak (2000) Frédéric Chyzak. 2000\.  An extension of Zeilberger’s fast algorithm to general holonomic functions.  _Discrete Math._ 217, 1-3 (2000), 115–134.  https://doi.org/10.1016/S0012-365X(99)00259-9
  * Cooley and Tukey (1965) James W. Cooley and John W. Tukey. 1965.  An algorithm for the machine calculation of complex Fourier series.  _Math. Comp._ 19 (1965), 297–301.  https://doi.org/10.2307/2003354
  * Couty et al. (2011) Danielle Couty, Jean Esterle, and Rachid Zarouf. 2011\.  Décomposition effective de Jordan-Chevalley.  _Gaz. Math._ 129 (2011), 29–49.
  * Czirbusz (2012) Sándor Czirbusz. 2012\.  Comparing the computation of Chebyshev polynomials in computer algebra systems.  _Ann. Univ. Sci. Budapest. Sect. Comput._ 36 (2012), 23–39.
  * Elaydi and Harris (1998) Saber N. Elaydi and William A. Harris, Jr. 1998.  On the computation of $A^{n}$.  _SIAM Rev._ 40, 4 (1998), 965–971.  https://doi.org/10.1137/S0036144597319235
  * Fiduccia (1985) Charles M. Fiduccia. 1985\.  An efficient formula for linear recurrences.  _SIAM J. Comput._ 14, 1 (1985), 106–112.  https://doi.org/10.1137/0214007
  * Flajolet and Salvy (1997) Philippe Flajolet and Bruno Salvy. 1997.  The SIGSAM challenges: symbolic asymptotics in practice.  _ACM SIGSAM Bull._ 31, 4 (1997), 36–47.  https://doi.org/10.1145/274888.274890
  * Gould (1972) Henry W. Gould. 1972\.  _Combinatorial identities_.  viii+106 pages.  A standardized set of tables listing 500 binomial coefficient summations.
  * Harvey and van der Hoeven (2021) David Harvey and Joris van der Hoeven. 2021.  Integer multiplication in time $O(n\log n)$.  _Ann. of Math. (2)_ 193, 2 (2021), 563–617.  https://doi.org/10.4007/annals.2021.193.2.4
  * Harvey and van der Hoeven (2022) David Harvey and Joris van der Hoeven. 2022.  Polynomial multiplication over finite fields in time $O(n\log n)$.  _J. ACM_ 69, 2 (2022), Art. 12, 40.  https://doi.org/10.1145/3505584
  * Hsieh et al. (1996) Po-Fang Hsieh, Mitsuhiko Kohno, and Yasutaka Sibuya. 1996\.  Construction of a fundamental matrix solution at a singular point of the first kind by means of the $SN$ decomposition of matrices.  _Linear Algebra Appl._ 239 (1996), 29–76.  https://doi.org/10.1016/S0024-3795(96)90003-8
  * Hyun et al. (2019) Seung Gyu Hyun, Vincent Neiger, and Éric Schost. 2019\.  Implementations of Efficient Univariate Polynomial Matrix Algorithms and Application to Bivariate Resultants. In _ISSAC’19_. ACM, 235–242.  https://doi.org/10.1145/3326229.3326272
  * Kauers and Mezzarobba (2019) Manuel Kauers and Marc Mezzarobba. 2019.  Multivariate Ore polynomials in SageMath.  _ACM Commun. Comput. Algebra_ 53, 2 (2019), 57–60.  https://doi.org/10.1145/3371991.3371998
  * Kauers and Pogudin (2017) Manuel Kauers and Gleb Pogudin. 2017.  Bounds for substituting algebraic functions into D-finite functions. In _ISSAC’17_. ACM, 245–252.  https://doi.org/10.1145/3087604.3087616
  * Koepf (1999) Wolfram Koepf. 1999\.  Efficient computation of Chebyshev polynomials in computer algebra.  In _Computer Algebra Systems: A Practical Guide_. Wiley, 79–99.
  * Koutschan (2010) Christoph Koutschan. 2010\.  A fast approach to creative telescoping.  _Math. Comput. Sci._ 4, 2-3 (2010), 259–266.  https://doi.org/10.1007/s11786-010-0055-0
  * Lipshitz (1988) L. Lipshitz. 1988\.  The diagonal of a $D$-finite power series is $D$-finite.  _J. Algebra_ 113, 2 (1988), 373–378.  https://doi.org/10.1016/0021-8693(88)90166-4
  * Massazza and Radicioni (2005) P. Massazza and R. Radicioni. 2005.  On computing the coefficients of bivariate holonomic formal series.  _Theoret. Comput. Sci._ 346, 2-3 (2005), 418–438.  https://doi.org/10.1016/j.tcs.2005.08.011
  * Miller and Brown (1966) J. C. P. Miller and D. J. Spencer Brown. 1966.  An algorithm for evaluation of remote terms in a linear recurrence sequence.  _Comput. J._ 9 (1966), 188–190.  https://doi.org/10.1093/comjnl/9.2.188
  * Schmidt (2000) Dieter Schmidt. 2000\.  Construction of the Jordan decomposition by means of Newton’s method.  _Linear Algebra Appl._ 314, 1-3 (2000), 75–89.  https://doi.org/10.1016/S0024-3795(00)00111-7
  * Shoup (2021) V. Shoup. 2021.  NTL: A library for doing number theory, v11.5.1.  https://libntl.org.
  * Stanley (1980) R. P. Stanley. 1980\.  Differentiably finite power series.  _European J. Combin._ 1, 2 (1980), 175–188.  https://doi.org/10.1016/S0195-6698(80)80051-5