Beating binary powering for polynomial matrices

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 complexity is \emph{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 \emph{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-finite sequence of polynomials, and modular exponentiation to the power $N$ for bivariate polynomials.

This can be achieved by binary powering for   , and in fact for   as well, since it is the top-right entry of   where  is the 2 × 2 companion matrix ( 0 1 1 1 ).This idea generalizes to any C-finite sequence (  ) ≥0 : a recurrence of order  ≥ 1 for (  ) ≥0 can be encoded, via its companion matrix, into an  × matrix recurrence of order 1.Then the term   of the sequence appears as the first entry of the product of the vector of initial values ( 0 , . . .,   −1 ) by the  th power of this  ×  companion matrix [18,31].Then   can be computed in  (log( )) arithmetic operations, and in  ( log( )) bit operations if (  ) ≥0 is an integer sequence, using fast integer multiplication [21].Here  is considered constant, i.e.,  ∈  (1).
Given  ∈ N, the direct iterative algorithm for computing   () has complexity  ( 2 ).It computes, for each  ≤  , all the  coefficients of the intermediate polynomial   (); in total this amounts to Θ( 2 ) coefficients.Therefore, if one wants to compute all of ( 0 , . . .,   ) 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   () alone.
To compute the polynomial   () faster, one can use, as in the scalar case, the reformulation of the second-order recurrence (1) as a first-order (polynomial) matrix recurrence: This shows that   () is the top-right entry of the matrix  ()  , where  () is the 2 × 2 companion matrix  () = ( 0 1 1  ).One can again compute  ()  using binary powering, whose costliest step is the multiplication of two polynomial matrices of degree about  /2.This yields   () in complexity  (M( )), where M( ) denotes the cost of polynomial multiplication in degree at most  .
Using FFT-based polynomial multiplication [11], this amounts to a number of operations in the base field K which is quasi-linear in  .Not only does this compare favorably to the complexity  ( 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 Θ( ) of the output polynomial   ().
In this context, the idea also generalizes to any C-finite sequence (  ()) ≥0 of polynomials in K[], which we will call polynomial C-finite sequences.Indeed, one can encode any recurrence of arbitrary (but independent of ) order  ≥ 1 and coefficients in K[] into a polynomial  ×  matrix recurrence of order 1, and the  th term of the sequence,   (), can be computed as an element in the  th power of an  ×  polynomial matrix multiplied by the polynomial vector of initial values.Conversely, computing the  th power of any polynomial matrix can be reduced to computing terms in polynomial C-finite sequences (see the introduction of Section 4).Binary powering allows to solve both problems in  (M( )) operations, and in  ( 2 log( )) bit operations if K = Q, considering both the recurrence order (or the matrix size)  and the recurrence degree (or the matrix degree)  as constant parameters, i.e., ,  ∈  (1).The main question addressed in this article is: Can one achieve a better complexity for these tasks?As far as scalar C-finite sequences are concerned, the algebraic complexity  (log( )) seems very difficult (if not impossible) to beat, but it is perhaps not impossible to improve the bit complexity  ( log( )) towards  ( ).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 Section 2. There, using Fibonacci polynomials as a test bench, we argue why it is indeed legitimate to hope for algorithms of complexity  ( ) for computing the  th term of a polynomial C-finite sequence.
Main result.Recall that a C-finite sequence is a sequence (  ) ≥0 of elements   in some ring  which satisfies a recurrence equation for  0 , . . .,   −1 ∈ .In this work we consider polynomial C-finite sequences, i.e., the case  = K[] for some (effective) field K of characteristic zero; thus   =   () ∈ K[].The customary data structure for representing such a sequence consists of the polynomials  0 (), . . .,   −1 () defining the recurrence and the  initial conditions  0 (), . . .,   −1 () ∈ K[].The order of the recurrence (3) is  while its degree is the maximum of the degrees of the   's.
Theorem 1.1.Let K be an effective field of characteristic 0. Let  and  be fixed positive integers.For each of the following problems, there exists an algorithm solving it in  ( ) operations (±, ×, ÷) in K: SeqTerm: Given a polynomial C-finite sequence (  ()) ≥0 of order and degree at most  and , compute the  th term   ().BivModPow: Given polynomials  (, ) and  (, ) in K[, ] of degrees in  and  at most  and , with  (, ) monic in , compute  (, )  mod  (, ).PolMatPow: Given a square polynomial matrix  () over K [𝑥] of size and degree at most  and , compute  ()  .
Our algorithms for these problems make essential use of divisions in K.We do not know if the complexity  ( ) can be achieved using only the operations (+, −, ×) in K.

Previous work.
As already mentioned, the classical way of computing the  th term of a given C-finite sequence uses binary powering of the companion matrix, see e.g.[31].Fiduccia's algorithm [18] utilizes binary powering in a polynomial quotient ring and improves the complexity with respect to  (but not with respect to  ).The fastest known algorithm [9] beats Fiduccia's by a constant factor.In the polynomial C-finite case and assuming ,  ∈  (1), all these algorithms have a complexity in  (M( )).
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   () = (−)   +1 (2) (with   () the th Fibonacci polynomial and  the imaginary unit) was considered in [27] (and later in [16]).These references present various methods for the computation of the Chebyshev polynomials (of the first and second kind) with complexity ranging from  ( ) to  ( 3 ).The results in [16,27] 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 [16,27] do not admit obvious generalizations.
An idea closely connected to a fundamental building block of our algorithms is explained in [19,Pbm. 4].There, Flajolet and Salvy exploit the fact that, given a polynomial  () in K[], the coefficient sequence of the th power  ()  satisfies a linear recurrence of order independent of , and with coefficients in K[, ] of degree independent of ; this recurrence allows them to compute (a selected coefficient of)  ()  more efficiently than by binary powering.This idea has been applied in [7, §8] to count points on hyperelliptic curves over finite fields, with applications to cryptography.The technique also yields a general solution to SeqTerm when  = 1.
Outline.The following observation generalizes that in [19]: the coefficient sequence of the th power of any algebraic function satisfies a recurrence of order and degree independent of .From this, in Section 3, we give algorithms for SeqTerm with cost  ( ).
To complete the proof of Theorem 1.1, we design reductions between the three problems.Obviously PolMatPow ⇒ SeqTerm, i.e., any algorithm for PolMatPow with cost  ( ) induces one for SeqTerm with cost  ( ) as well.Indeed, the  th term of a polynomial C-finite sequence is equal to an entry of the product of the vector of initial values and the  th power of a companion matrix, and this polynomial vector-matrix multiplication costs  ( ).Conversely, it also holds that SeqTerm ⇒ PolMatPow.One natural way to see this is to consider  2 sequences corresponding to each entry of  ()  , with recurrence given by the characteristic polynomial of  (); see the introduction of Section 4. In Section 4.2, we give a more efficient algorithm for this reduction, based on an algorithm for SeqTerm ⇒ BivModPow described in Section 4.1.
Basics of complexity and D-finite functions.Hereafter, K denotes an effective field of characteristic zero.We analyze the performance of algorithms in the algebraic complexity model, meaning that arithmetic operations (±, ×, ÷) in the base field K are counted at unit cost.As before, M( ) stands for the complexity of multiplying two polynomials in K[] of degree at most  .With FFT-based multiplication M( ) ∈  ( log( ) log log( )) [11], improved to  ( log( )) if K contains suitable roots of unity [14] or if K is a finite field [22].A power series  () ∈ K[[]] is said to be D-finite if it satisfies a linear differential equation (LDE) of the form for some  0 (), . .

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  th Fibonacci polynomial   () in complexity  ( ).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  (, ) ] of the sequence (  ()) ≥0 is rational, and equal to /(1 −  −  2 ).
As mentioned in the introduction, the analogue of formula (8) for the case of Chebyshev polynomials of the first kind   () was already exploited in [27, §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).

Second method via algebraic substitution
There is another method for computing   () in  ( ), which has the advantage that it generalizes to any C-finite sequence, as we will show in Section 3.1.The crucial remark (Lemma 3.2) is that since  ± () is algebraic,  ± ()  satisfies a "small" LDE, of order and degree independent of .The same holds for 1/( + () − − ()), therefore for   () as well.More precisely,  ± ()  satisfies the LDE Using (5), it then follows that the polynomial   () satisfies Writing   () =  −1 =0     , plugging into (10) for  =  and extracting the ( + 2)nd coefficient, it now follows that the sequence (  )  ≥0 satisfies recurrence (9).The initial conditions  0 ,  1 are given by   () mod  2 which can be found in complexity  (log( )) by computing the  th power of the companion matrix (2) in K[]/( 2 ) by binary powering and reducing mod  2 in each step.As before, unrolling recurrence ( 9) with these initial terms provides a way to compute   () in complexity  ( ).

Third method via Creative Telescoping
Writing  (, ) = /(1 −  −  2 ) we are interested in a differential equation for the coefficient of   in  (, ).By Cauchy's integral formula, we have for sufficiently small  > 0: Then the method of creative telescoping [1] can be used to find an LDE 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 . By Cauchy's integral theorem, the contour integral of the right-hand side vanishes, and ( 10) follows.Then one can conclude in the same way as in the previous method and compute   () in complexity  ( ).

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   () in complexity  ( ) is the existence of the recurrence (9) (or equivalently the LDE (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   () is a polynomial, it does satisfy the first-order LDE   () ′ () −  ′  () () = 0, but this one is trivial for our purposes.Indeed, converting this differential equation into a recurrence satisfied by the sequence of coefficients of   () yields a recurrence of order deg(  ), which is obviously useless for computing the coefficients of   .Rather, we would like to find an LRE/LDE whose order and degree are independent of  .This is the purpose of the next section.Specifically, in §3.1 we explain how it can be computed by algebraic substitution (generalizing §2.2) and in §3.2 we show that it can also be found via creative telescoping (generalizing §2.3).

Computing 𝑢 𝑁 (𝑥) in 𝑂 (𝑁 )
By generalizing the ideas of Section 2.2, it is not difficult to prove that the th term of a polynomial C-finite sequence (  ()) ≥0 satisfies an LDE whose order and degree are independent of , and consequently, that there exists a linear recurrence relation for the coefficient sequence of   () whose order (say ) and degree are again independent of .Then, for a given  ∈ N, first computing initial terms by binary powering of the companion matrix in K[]/(  ) and then unrolling this recurrence for  =  , we achieve a complexity  ( ) for the computation of   ().
Theorem 3.1.Let (  ()) ≥0 be a polynomial C-finite sequence.Then there exists   ∈ K[, ]⟨  ⟩ with order and degree independent of , and such that   (  ()) = 0. Consequently, writing   () =  ≥0  ,   , there exist, for some  ∈ N independent of , polynomials  0 (, ), . . .,   (, ) ∈ K[, ] of degrees independent of , and such that the sequence ( , )  ≥0 satisfies the recurrence In the theorem above it is crucial that neither the order nor the degree of   depend on .Since each   () is a polynomial, it is a tautology to say that it satisfies some LDE: one may simply take  =    , where  > deg(  ()) or  =   ()  −  ′  ().However, it is a nontrivial fact that   () satisfies an LDE of the form The most direct proof of Theorem 3.1 uses the explicit expression (13) for   () and the following classical fact about algebraic substitution into D-finite functions.Recall that a function () is called algebraic over K() if it satisfies a nontrivial polynomial relation  (, ()) = 0 for some  (, ) ∈ K[, ].Size and complexity bounds on differential equations for algebraic functions, and more generally on algebraic substitution, are given in [5,26].Lemma 3.2.Let () be an algebraic function over K() and let () be D-finite.Then  () = (()) is D-finite.In particular, ()  satisfies an LDE of order and degree independent of .
Proof.The first part is a classical result, see for example [34,Thm. 2.7].In the proof one shows that the vector space spanned over K() by ( ( ) Proof of Theorem 3.1.Write   () as in (13).By Lemma 3.2, each   ()  satisfies an LDE of order and degree independent of , hence the same holds for   (, )  ()  , and finally for   ().It follows that the coefficient sequence of   () is P-finite with order and degree independent of .□ Since all steps in the proofs above are effective and independent of  , this leads to Algorithm 1.Its Lines 2 to 7 can be seen as "precomputations" since they do not depend on  .As already mentioned, Line 8 has complexity  (log( )) and Line 9 has complexity  ( ).Thus, Algorithm 1 solves SeqTerm in complexity  ( ), up to a potential issue during the unrolling at Line 9 of the recurrence from Line 7. Indeed, this unrolling may be impossible for some values , namely those for which   ( , ) vanishes.We will explain how to overcome this problem in Section 3.3.
For practical applications, however, computing the polynomials   (, ) in Line 4 as well as the LCLM in Line 6 is algorithmically somewhat cumbersome.Thus, generalizing the approach in Section 2.3, we now propose a variant of Algorithm 1 which replaces Lines 2 to 6 by an algorithm based on creative telescoping.
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 [2].Let L = K().=0   ()   is a telescoper for  .This procedure cannot be directly applied to  (, )/ +1 if  is an indeterminate.At the same time, if  =  ∈ N is fixed, it is a priori not obvious that deg    () will be independent of  .Moreover, the complexity of the algorithm will depend on  , which we want to avoid.As we will now explain, to achieve this, one should see  (, )/ +1 not as a rational function in  and  with potentially large degree in the numerator, but as a hyperexponential function with the parameter  appearing solely as a coefficient in the logarithmic derivative.Recall that  (, ) is called hyperexponential if both    / and    / belong to K(, ).
For hyperexponential functions, the Almkvist-Zeilberger algorithm [1] 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 th Fibonacci polynomial for a variable .Note that if  is specialized to an integer  before the execution of the command above, the implemented algorithm becomes slower as  grows.
It is, however, not clear that the Almkvist-Zeilberger algorithm applied to  (, )/ +1 will always find a telescoper whose degree and order are independent of , even though we know from Section 3.1 that an LDE 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 [3].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 LDE for   ().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  (, )/ +1 , and consequently an LDE for   (), whose order and degree do not depend on .
where  1 (, ) is  1 () with the th coefficient   replaced by   /( − ).Setting  =  1 and  =  1 proves the induction basis.Now assume that  − > 0 and note that , so equation ( 16) is equivalent to The Hermite reduction applied to Comparing ( 17) and ( 18 the first integral vanishes by Cauchy's integral theorem, and the second integral vanishes by construction of the   (, ).This provides a variant for Lines 2 to 6 of Algorithm 1, as described in Algorithm 3. The above-mentioned potential issue with unrolling persists; the next section deals with this problem.

The singular case
In this section we discuss the potential issue of our algorithm that can occur if the sequence for the coefficients of   () cannot be unrolled due to singularities.We shall first highlight this problem and its solution by means of an example, then in the last paragraph of this section we explain the general strategy.
Consider the polynomial C-finite sequence   () given by for all  ≥ 0 with initial conditions  0 = 3,  1 = ).With the initial conditions and after a partial fraction decomposition it follows that the generating function of   () is given by Hence, the solution is   () = 2  +   +  2 and can be written down in  (log( )) operations.However, as we shall explain now, the direct application of any of the methods described earlier fails.
According to Theorem 3.1,   () satisfies an LDE whose degree and order are independent of .Indeed, using creative telescoping one quickly finds an annihilator for   () = ∮  (, )/ +1 d: In other words,  , = 0 for all  ∈ N except  ∈ {0, , 2}.In order to "unroll" this recurrence we need to know  ,0 ,  , and  ,2 .However, it is not immediately clear how to compute those terms for  =  in  ( ) operations from the initial input (without using the explicit solution).
We propose the following easily generalizable solution: consider   () =   ( + 1).Then the LDE for   () is given by ( + 1) 2  3  − 3( + 1)( − 1) 2  + (2 − 1)( − 1)    () = 0, and for the coefficient sequence of   () =  ≥0  ,   we find Now the leading coefficient of the recurrence is ( + 1) ( + 2) ≠ 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  2 .Having computed   (), it remains to find   () =   ( − 1).Note that expanding the polynomial results in an  (M( )) algorithm.However, recall from ( 19) that we only need to compute   , and   ,2 , or, in other words, the coefficients of   and  2 in   ( − 1).For any  it holds that and the sum is finite because   () is a polynomial.Clearly, it can be computed in complexity  ( ) for any .
Generally speaking, an issue with unrolling the recurrence for ( , )  ≥0 occurs if the roots of the leading polynomial are positive integers that depend on .Indeed, roots that are nonintegral clearly do not cause any problems in the unrolling step and if a root is independent of  then we may just compute more initial terms while the complexity of this step stays bounded by  (log  ).Let  be the set of the problematic roots.Note that the size of  is independent of  since  is a subset of the roots of the leading polynomial   (, ) in ( 14) and deg    (, ) is bounded independently of  by Theorem 3.1.Moreover,  can be nonempty only if the LDE for   () is singular at 0 (that is, if  ℓ () =  ℓ (, ) in ( 4) vanishes at  = 0).In this case, one can always define   () =   ( + ) for  ∈ K a nonsingular point of the LDE ( ℓ (, ) ≠ 0).Then the coefficients  , of   () can be computed from  (1) initial conditions via unrolling a recurrence.Using the formula (20) (with − instead −1) and the fact that   () is a polynomial, one can compute the coefficients   , for  ∈ .It is then possible to unroll the recurrence for (  , )  ≥0 and find   () in complexity  ( ).

Computing bivariate modular powers
Let L = K[] and ,  ∈ L[].Assume that , seen as a univariate polynomial in  of degree  , is monic.For  ∈ N, the Euclidean division in L[] ensures the existence of unique ,  ∈ L[] such that deg  () <  and   =  +.The polynomial  is   mod .Assume that  and  are fixed, and let  deg  () (which is thus in  (1)).Then, writing  =  −1 =0   ()  , it holds that deg    () =  ( ).The efficient computation of  when  = , given  (, ) and  , is the first step for proving BivModPow in Theorem 1.1.
This strategy, outlined on an example, generalizes in the obvious way.Explicitly, we have the following lemma (see [9,Lem. 2]).

EXPERIMENTS
The main precomputation step for all our algorithms consists in starting with a rational function  ∈ K(, ) ∩ K[[, ]] and in finding a differential operator   that annihilates   () = [  ] (, ) and whose degree and order are independent of .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  (, )/ +1 is the reductionbased creative telescoping in Maple.More specifically, redct is the fastest, followed by HT.The implementation in ore_algebra [25] 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,  and  are large and  is small; this is expected since these ranges correspond to cases where the order of the recurrence to be unrolled is close to  .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 to 3; as expected it grows when  grows, with  and  fixed.
For large  , 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  = 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.

PERSPECTIVES
We have shown that it is possible to beat, both in theory and in practice, the 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, left for future investigations.
More detailed complexity analysis.The 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  and .For simplicity, these parameters were assumed to be  (1) in this work.For the  th power of an  × matrix  () of degree , binary powering has complexity  (M() 2 +   ), where  ∈ [2, 3] is a feasible exponent of matrix multiplication over K.With our approach, it is legitimate to target a differential equation satisfied by the entries of  ()  of order  with coefficients in  of degree  ( 3 ), yielding a recurrence of order  ( 3 ) and coefficients in  of degree at most  .For large  , this would result in a complexity in  ( 2 M( )).Using different LDEs, of order  ( ) and coefficients of degree  ( 2 ) could even lead to  ( M( )).
The th coefficient of the  th term.For some (large) integers  ,  ∈ N, one might be interested in computing the single coefficient [    ] (, ) of a rational function Equivalently it is natural to wonder: how fast can one compute the th coefficient of the  th term of a C-finite sequence (  ()) ≥0 ?Using our method, a recurrence with initial conditions for the coefficients of   () can be deduced in  (log( )) operations.Then (assuming that the recurrence is nonsingular) the th coefficient can be found in  (M( √ )) operations by using babysteps/giant-steps techniques [7,12].We expect that, at least under a genericity assumption, this problem can be solved in complexity  (log( ) + M( √ )) which is a big improvement compared to the previous best  ( + ) by Massazza and Radicioni [30].
Polynomial P-finite sequences.A somewhat related task is to study the analogous problem to SeqTerm for polynomial P-finite sequences, that is for We expect that, at least under a genericity assumption, a generalization of Lemma 3.4 (based on results in [4,35]) should exist, implying in particular that   () satisfies an LDE of order and degree independent of  .Generalizing this even further, one might study the creative telescoping problem for rational functions of the form  (x) =  ( 1 ,...,  )  ( 1 ,...,  ) ( 1 ,...,  )  .We expect that (at least generically) the minimal telescoper for  (x) has order and degree independent of  and can be found via a Griffiths-Dwork reduction type approach, based on ideas from [8].
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.[15,17,23,32].It ensures that any polynomial matrix  ∈ K[]  × can be written as  =  + , where  ∈ K()  × is diagonalizable over K(),  ∈ K()  × is nilpotent, and  = .From this decomposition it follows that   =  −1 =0     −   .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  ( ) 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  (, ) = /(1 −  −  2 ) of   () satisfies the linear PDE and extracting the coefficient of     in (22) [29,Lem. 3] proves that such a PDE exists for any rational function  (, ).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 [1] or reduction-based [2]), and how the corresponding method compares to the aforementioned ones.
Integer case in bit complexity  ( ).Recall the analogy between the bit complexity for finding the  th term of a C-finite sequence over Z and the complexity for finding the  th term of a C-finite sequence over K[].Our work achieves  ( ) for the latter, so it is now natural to target  ( ) for the former, for instance for the  th Fibonacci number or simply 3  .This remains widely open.