Alin Bostan
email: alin.bostan@inria.fr

Ryuhei Mori
email: mori@c.titech.ac.jp

A Simple and Fast Algorithm for Computing the N -th Term of a Linearly Recurrent Sequence

Keywords: Algebraic Algorithms, Computational Complexity, Linearly Recurrent Sequence, Rational Power Series, Fast Fourier Transform

, has the same arithmetic complexity up to a constant factor. Our algorithm is simpler, faster and obtained by a totally different method. We also discuss several algorithmic applications, notably to polynomial modular exponentiation and powering of matrices.

Introduction

1.1 General context Computing efficiently selected terms in sequences is a basic and fundamental algorithmic problem, whose applications are ubiquitous, for instance in theoretical computer science [START_REF] Ravikumar | Weak minimization of DFA-an algorithm and applications[END_REF][START_REF] Lürwer-Brüggemeier | On faster integer calculations using non-arithmetic primitives[END_REF], algebraic complexity theory [START_REF] Paterson | On the number of nonscalar multiplications necessary to evaluate polynomials[END_REF][START_REF] Strassen | Polynomials with rational coefficients which are hard to compute[END_REF], computer algebra [START_REF] Giesbrecht | Nearly optimal algorithms for canonical matrix forms[END_REF][START_REF] Storjohann | High-order lifting[END_REF][START_REF] Lecerf | New recombination algorithms for bivariate polynomial factorization based on Hensel lifting[END_REF], cryptography [START_REF] Gong | Public-key cryptosystems based on cubic finite field extensions[END_REF][START_REF] Gong | The GH publickey cryptosystem[END_REF][START_REF] Giuliani | New LFSR-Based Cryptosystems and the Trace Discrete Log Problem (Trace-DLP)[END_REF][START_REF] González | Discrete logarithm like problems and linear recurring sequences[END_REF], algorithmic number theory [START_REF] Strassen | Einige Resultate über Berechnungskomplexität[END_REF][START_REF] Agrawal | PRIMES is in P[END_REF], effective algebraic geometry [START_REF] Bostan | Linear recurrences with polynomial coefficients and application to integer factorization and Cartier-Manin operator[END_REF][START_REF] Harvey | Counting points on hyperelliptic curves in average polynomial time[END_REF], numerical analysis [START_REF] Mezzarobba | NumGfun: a package for numerical and analytic computation and D-finite functions[END_REF][START_REF] Mcnamee | Numerical methods for roots of polynomials[END_REF] and computational biology [START_REF] Nuel | Sparse approaches for the exact distribution of patterns in long state sequences generated by a Markov source[END_REF].

In simple terms, the problem can be formulated as follows:

Given a sequence (u n ) n≥0 in an effective ring 0 R, and given a positive integer N ∈ N, compute the term u N as fast as possible.

Here, the input (u n ) n≥0 ∈ R N is assumed to be a recurrent sequence, specified by a data structure 0 The ring R is assumed to be commutative with unity and effective in the sense that its elements are represented using some data structure, and there exist algorithms for performing the basic ring operations (+, -, ×) and for testing equality of elements in R.

consisting in a recurrence relation and sufficiently many initial terms that uniquely determine its terms.

Efficiency is measured in terms of ring operations (algebraic model), or of bit operations (Turing machine model). The cost of an algorithm is respectively estimated in terms of arithmetic complexity or of binary complexity. Both measures have their own usefulness: the algebraic model is relevant when ring operations have essentially unit cost (typically, if R is a finite ring such as the prime field F p := Z/pZ), while the bit complexity model is relevant when elements of R have a variable bitsize, and thus ring operations in R have variable cost (typically, when R is the ring Z of integer numbers).

The recurrence relation satisfied by the input sequence (u n ) n≥0 might be of several types: (C) linear with constant coefficients, that is of the form,

u n+d = c d-1 u n+d-1 + • • • + c 0 u n , n ≥ 0,
for some given c 0 , . . . , c d-1 in R. In this case we simply say that (u n ) n≥0 is linearly recurrent (or, C-recursive). The most basic examples are the geometric sequence (q n ) n≥0 , for q ∈ R, and the Fibonacci sequence (F n ) n with F n+2 = F n+1 + F n for n ≥ 0 and F 0 = 0, F 1 = 1.

(P) linear with polynomial coefficients, of the form,

u n+d = c d-1 (n)u n+d-1 + • • • + c 0 (n)u n , n ≥ 0,
for some given rational functions c 0 (x), . . . , c d-1 (x) in R(x). In this case the sequence is called holonomic (or, P-recursive). Among the most basic examples, other than the C-recursive ones, there is the factorial sequence (n!) n≥0 = (1, 1, 2, 6, 24, 120, . . .) and the Motzkin sequence (u n ) n≥0 = (1, 1, 2, 4, 9, 21, 51, . . .) specified by the recurrence u n+1 = 2n+3 n+3 • u n + 3n n+3 • u n-1 and the initial conditions u 0 = u 1 = 1.

(Q) linear with polynomial coefficients in q and q n , that is of the form,

u n+d = c d-1 (q, q n )u n+d-1 +• • •+c 0 (q, q n )u n , n ≥ 0,
for some q ∈ R and some rational functions c 0 (x, y), . . . , c d-1 (x, y) in R(x, y). In this case, the sequence is called q-holonomic; such a sequence can be seen as a q-deformation of a holonomic sequence (in the sense that when q → 1, the limit sequence tends to be holonomic). A typical example is the q-factorial [n] q ! := (1 + q)• • •(1 + q + • • • + q n-1 ).

In all these classes of examples, the recurrence is linear, and the integer d that governs the length of the recurrence relation is called the order of the corresponding linear recurrence.

Of course, some interesting sequences satisfy nonlinear recurrences, as is the case for the so-called Somos-4 sequence (1, 1, 1, 1, 2, 3, 7, 23, . . .) defined by: u n+4 = u n+3 u n+1 +u 2 n+2 /u n together with u 0 = • • • = u 3 = 1, but we will not consider this larger class in what follows.

To compute the N -th term u N of a sequence of type (P), resp. (Q), the best known algorithms are presented in [START_REF] Chudnovsky | Approximations and complex multiplication according to Ramanujan[END_REF][START_REF] Bostan | Linear recurrences with polynomial coefficients and application to integer factorization and Cartier-Manin operator[END_REF], resp. in [START_REF] Bostan | Computing the N -th Term of a q-Holonomic Sequence[END_REF]. In the algebraic model, they rely on an algorithmic technique called baby-step/giant-step, allowing to compute u N in a number of operations in R that is almost linear in √ N , up to logarithmic factors. This should be contrasted with the direct iterative algorithm, of arithmetic complexity linear in N .

In the bit model, the same references provide different algorithms based on a different technique, called binary splitting; these algorithms are quasi-optimal in the sense that they are able to compute u N in a number of bit operations almost linear (up to logarithmic factors) in the bitsize of the output value u N . Once again, this should be contrasted with the direct iterative algorithm, whose binary complexity is larger by at least one order of magnitude (e.g., in case (P) the naive algorithm has bit complexity O(N 3 )).

The case of C-recursive sequences

In what follows, we will restrict our attention to the case (C) only. This case obviously is a subcase of both cases (P) and (Q). It presents an exceptional feature with respect to the algebraic model: contrary to the general cases (P) and (Q), in case (C) it is possible to compute the term u N using a number of arithmetic operations in R that is only logarithmic in N .

For the geometric sequence u n = q n , this is known since Pingala (∼200 BC) who seemingly is the inventor of the algorithmic method of binary powering, or squareand-multiply [46, §4.6.3]. The corresponding algorithm is recursive and based on the equalities

q N = (q N/2 ) 2 , if N is even, q • (q N -1 2 ) 2 , else.
The arithmetic complexity of this algorithm is bounded by 2 log N multiplications 1 in R, which represents a tremendous improvement compared to the naive iterative algorithm that computes the term q N in N -1 multiplications in R, by simply unrolling the recurrence u n+1 = q • u n with q 0 = 1.

In the general case (C), Miller and Spencer Brown showed in 1966 [START_REF] Miller | An algorithm for evaluation of remote terms in a linear recurrence sequence[END_REF] that a similar complexity can be obtained by converting 2 the scalar recurrence of order d

(1.1) u n+d = c d-1 u n+d-1 + • • • + c 0 u n , n ≥ 0, into a vector recurrence of order 1 (1.2)      u n u n+1 . . . u n+d-1      vn =      1 . . . 1 c 0 c 1 • • • c d-1      M ×      u n-1 u n . . . u n+d-2      vn-1 , n ≥ 1,
and by using binary powering in the ring M d (R), of d×d matrices with coefficients in R, in order to compute M N recursively by

M N = (M N/2 ) 2 , if N is even, M • (M N -1 2 ) 2 , else.
From there, u N can be read off the matrix-vector product v N = M N • v 0 . The complexity of this method is O(d θ log N ) operations in R, where θ ∈ [2, 3] is any feasible exponent for matrix multiplication in M d (R).

Strangely enough, the paper [START_REF] Miller | An algorithm for evaluation of remote terms in a linear recurrence sequence[END_REF] of Miller and Spencer Brown was largely overlooked in the subsequent literature, and their result has been rediscovered several times in the 1980s. For instance, Shortt [START_REF] Shortt | An iterative program to calculate Fibonacci numbers in O(log n) arithmetic operations[END_REF] proposed a O(log N ) algorithm for computing the N -th Fibonacci number 3,4 and extended it together with Wil- 1 In this article, notation log refers to the logarithm in base 2. 2 [15, p. 74] calls this un truc bien connu ("a well-known trick"). 3 Shortt's algorithm had actually appeared before, in the 1969 edition of Knuth's book [45, p. 552], as a solution of Ex. 26 (p. 421, §4.6.3). The algorithm is based on the so-called doubling formulas

(F 2n , F 2n-1 ) = (F 2 n +2FnF n-1 , F 2 n +F 2 n-1
), actually due to Lucas (1876) and Catalan (1886), see e.g. [START_REF] Dickson | History of the theory of numbers[END_REF]Ch. XVII]. The currently best implementation for computing F N over Z (mpz_fib_ui from GMP) uses a variant of this method, requiring just two squares (and a few additions) per binary digit of N . 4 Already in 1899, G. de Rocquigny asked "for an expeditious procedure to compute a very distant term of the Fibonacci sequence" [START_REF] De Rocquigny | Question 1541. [I25a[END_REF]. In response, several methods (including the one mentioned by Knuth in [45, p. 552]) have been published one year later by Rosace (alias), E.-B. Escott, E. Malo, C.-A. Laisant and G. Picou [START_REF] Rosace | Answers to Question 1541[END_REF]. This fact does not seem to have been noticed in the modern algorithmic literature before the current paper, although the reference [START_REF] Rosace | Answers to Question 1541[END_REF] is mentioned in Dickson's formidable encyclopedic book [20, p. 404].

son [START_REF] Wilson | An O(log n) algorithm for computing general order-k Fibonacci numbers[END_REF] to the computation of order-d Fibonacci numbers in O(d 3 log N ) arithmetic operations. The same cost has also been obtained by Dijkstra [START_REF] Dijkstra | In honour of Fibonacci[END_REF] and Urbanek [START_REF] Urbanek | An O(log n) algorithm for computing the nth element of the solution of a difference equation[END_REF]. Pettorossi [START_REF] Pettorossi | Derivation of an O(k 2 log n) algorithm for computing order-k Fibonacci numbers from the O(k 3 log n) matrix multiplication method[END_REF], and independently Gries and Levin [START_REF] Gries | Computing Fibonacci numbers (and similarly defined functions) in log time[END_REF], improved the algorithm and lowered the cost to O(d 2 log N ), essentially by taking into account the sparse structure of the matrix M . See also [START_REF] Er | Computing sums of order-k Fibonacci numbers in log time[END_REF][START_REF] Martin | A presentation of the Fibonacci algorithm[END_REF][START_REF] Er | A formal derivation of an O(log n) algorithm for computing Fibonacci numbers[END_REF][START_REF] Er | An O(k 2 log(n/k)) algorithm for computing generalized order-k Fibonacci numbers with linear space[END_REF][START_REF] Holloway | Algorithms for computing Fibonacci numbers quickly[END_REF][START_REF] Protasi | On the number of arithmetical operations for finding Fibonacci numbers[END_REF][START_REF] Giuliani | A new algorithm to compute remote terms in special types of characteristic sequences[END_REF][START_REF] Khomovsky | Efficient computation of terms of linear recurrence sequences of any order[END_REF] for similar algorithms.

Fiduccia's algorithm

The currently best algorithm is due to Fiduccia5 [START_REF] Fiduccia | An efficient formula for linear recurrences[END_REF]. It is based on the following observation: the matrix M in (1.2) is the transpose of the companion matrix C which represents the R-linear multiplication-by-x map from the quotient ring

R[x]/(Γ) to itself, where Γ = x d - d-1 i=0 c i x i .
Therefore, denoting by e the row vector e = [1 0

• • • 0], the N -th term u N equals (1.3) u N = e•v N = e•M N •v 0 = C N • e T T •v 0 = x N mod Γ, v 0 ,
where the inner product takes place between the vector

v 0 = u 0 • • • u d-1
of initial terms of (u n ) n≥0 , and the vector whose entries are the coefficients of the remainder (x N mod Γ) of the Euclidean division of x N by Γ. Thus, computing u N is reduced to computing the coefficients of (x N mod Γ), and this can be performed efficiently by using binary powering in the quotient ring In conclusion, Fiduccia's algorithm allows the computation of the N -th term u N of a linearly recurrent sequence of order d using O(M(d) log N ) operations in R. Since 1985, this is the state-of-the-art algorithm for this task in case (C).

A := R[x]/(Γ),
A closer inspection of the proof of [27, Corollary 9.7] shows that a more precise estimate for the arithmetic cost of Fiduccia's algorithm is

(1.4) F(N, d) = 3 M(d) log N + O(d log N )
operations in R. This comes from the fact that squaring

6 in A = R[x]/(Γ)
T(N, d) = 2 M(d) log(N + 1) + M(d)
arithmetic operations in R.

The proof of this result is based on a very natural algorithm, which will be presented in Section 2.1. Let us remark that it improves by a factor of 1.5 the complexity of Fiduccia's algorithm. This factor is even higher in the FFT setting, where polynomial multiplication is assumed to be performed using Fast Fourier Transform techniques. In this setting, we obtain the following complexity result, which will be proved in Section 4.

Theorem 2. One can compute the N -th term of a linearly recurrent sequence of order d with coefficients in a field K supporting FFT using7

∼ 2 3 M(d) log(N )
arithmetic operations in K.

Algorithms 1 and 2 (underlying Theorem 1) and Algorithm 11 (underlying Theorem 2) are both of LSBfirst (least significant bit first) type. This prevents them from computing simultaneously several consecutive terms of high indices, such as u N , . . . , u N +d-1 . This makes a notable difference with Fiduccia's algorithm from §1.3. For this reason, we will design a second algorithm, of MSB-first (most significant bit first) type, by "transposing" Algorithm 1.

This leads to the following complexity result.

Theorem 3. One can compute the terms of indices

N -d+1, . . . , N of a linearly recurrent sequence of order d with coefficients in a ring R using

(2 M(d) + d) log(N + 1) + O(M(d))
arithmetic operations in R.

The method underlying this complexity result is based on Algorithms 5, 6 and 8, which are presented in Section 3. Along the way, using the MSB-first Algorithm 5, we improve the cost of polynomial modular exponentiation, which is a central algorithmic task in computer algebra, with many applications. Since this result has an interest per se, we isolate it here as our last complexity result.

Theorem 4. Given N ∈ N and a polynomial Γ(x) in R[x] of degree d, one can compute x N mod Γ(x) using (2 M(d) + d) log(N + 1) + M(d) arithmetic operations in R.
This complexity of ∼ 2 M(d) log N compares favorably with the currently best estimate of ∼ 3 M(d) log N obtained by square-and-multiply in the quotient ring R[x]/(Γ(x)), combined with fast modular multiplications performed either classically [27, Corollary 9.7], or using Montgomery's algorithm [START_REF] Montgomery | Modular multiplication without trial division[END_REF]. In the FFT setting, the gain is even larger, and our results improve on the best estimates, due to Mihăilescu [START_REF] Mihăilescu | Fast convolutions meet Montgomery[END_REF].

Structure of the paper

In Section 2 we propose our LSB-first (least significant bit first) algorithm for computing the N -th term of a C-recursive sequence. We design in Section 3 a second algorithm, which is an MSB-first (most significant bit first) variant, and discuss several algorithmic applications, including polynomial modular exponentiation and powering of matrices. In Section 4 we specialize and analyze Algorithm 1 in the specific FFT setting, where polynomial multiplication is based on Discrete Fourier Transform techniques, and we compare it with the FFT-based Fiduccia's algorithm. We conclude in Section 5 by a summary of results and plans of future work.

The LSB-first algorithm and applications

We will prove Theorem 1 in §2.1, where we propose the first main algorithms (Algorithms 1 and 2), which are faster than Fiduccia's algorithm. Then, in §2.2 we instantiate them in the particular case of the Fibonacci sequence. The resulting algorithm is competitive with state-of-the-art algorithms.

LSB-first algorithm: Proof of Theorem 1

In this section, we give the proof of Theorem 1. Let us denote by

F (x) ∈ R[[x]] the generating function of the sequence (u n ) n≥0 , F (x) := n≥0 u n x n .
We use the following classical characterization of generating functions of linearly recurrent sequences [START_REF] Ranum | The general term of a recurring series[END_REF]. To be self-contained, we include a proof.

Lemma 1. A sequence (u n ) n≥0 is linearly recurrent of order d if and only if its generating function is

F (x) = P (x)/Q(x) for some polynomials P (x) of degree at most d -1 and Q(x) of degree d satisfying Q(0) = 1. Proof. Assume (u n ) n≥0 is a linearly recurrent sequence of order d. Let us define Q(x) := 1 -c d-1 x -• • • -c 0 x d , that is the reversal of the characteristic polynomial Γ(x) = x d - d-1 i=0 c i x i of recurrence (1.1). Then, P (x) := Q(x) • F (x)
is a polynomial of degree less than d. This is immediately seen by checking that, for any n ≥ 0, the coefficient of x n+d in the power series

Q(x) • F (x) is equal to u n+d -c d-1 u n+d-1 -• • • -c 0 u n , hence it is zero by (1.1).
Conversely, assume that the generating function of a sequence (u n ) n≥0 is F (x) = P (x)/Q(x) for some polynomials P (x) and Q(x) satisfying the conditions in the lemma. Then, Q(x)F (x) = P (x) implies that the coefficient of x n+d in Q(x)F (x) is zero for all integers n ≥ 0. Hence, the linear recurrence equation (1.1) must be satisfied, where 1

-c d-1 x -• • • -c 0 x d := Q(x). Proof of Theorem 1. From the values u 0 , . . . , u d-1 and c 0 , . . . c d-1 , we obtain Q(x) := 1 -c d-1 x -• • • -c 0 x d without any operations in R, and can compute P (x) := (u 0 + • • • + u d-1 x d-1 ) • Q(x) mod x d using at most M(d) operations in R. According to Lemma 1, our goal is to compute u N = [x N ] P (x) Q(x) , where [x N ] F (x) denotes the coefficient of x N in F (x) ∈ R[[x]]. From Lemma 1, we can assume that Q(0) = 1.
But, for the sake of generality, in the following, we only assume that x) by the numerator P (x) and the denominator Q(x), and obtain

Q(0) is invertible in R. If N = 0, we have u N = P (0)/Q(0). If N ≥ 1, we multiply Q(-
u N = [x N ] P (x)Q(-x) Q(x)Q(-x) .
Here, Q(x)Q(-x) is an even polynomial, which means that all coefficients of x 2n+1 for integers n ≥ 0 in this Algorithm 1 (OneCoeff)

Input: P (x), Q(x), N Output: [x N ] P (x) Q(x) Assumptions: Q(0) invertible and deg(P ) < deg(Q) =: d 1: while N ≥ 1 do 2: U (x) ← P (x)Q(-x) U = 2d-1 i=0 U i x i 3:
if N is even then

4: P (x) ← d-1 i=0 U 2i x i 5:
else 6:

P (x) ← d-1 i=0 U 2i+1 x i 7: A(x) ← Q(x)Q(-x) A = 2d i=0 A i x i 8: Q(x) ← d i=0 A 2i x i 9: N ← N/2 10: return P (0)/Q(0)
polynomial are equal to 0. Hence, there exists a unique polynomial

V (x) satisfying V (x 2 ) = Q(x)Q(-x).
Here, the degree of V (x) is d and

V (0) = Q(0) 2 is invertible in R. Now, we obtain u N = [x N ] U (x) V (x 2 )
for U (x) := P (x)Q(-x). The numerator U (x) has degree at most 2d -1. Here, 1/V (x 2 ) is an even formal power series. That implies that if N is even (or odd), we can ignore the odd (or even) part of U (x). Let U e (x) and U o (x) be the even and odd parts of U (x), respectively, i.e., U e (x) :=

d-1 i=0 U 2i x i and U o (x) := d-1 i=0 U 2i+1 x i for U (x) = 2d-1 i=0 U i x i .
Then, we obtain the decomposition U (x) = U e (x 2 ) + xU o (x 2 ). Hence,

u N = [x N ] Ue(x 2 ) V (x 2 ) , if N is even [x N ] xUo(x 2 ) V (x 2 ) , else. = [x N/2 ] Ue(x) V (x) , if N is even [x (N -1)/2 ] Uo(x) V (x) , else.
Here, U e (x) and U o (x) are polynomials of degree at most d -1, while V (x) is a polynomial of degree d satisfying that V (0) is invertible in R. Hence, we can repeat this reduction until N ≥ 1.

Our algorithm for computing [x N ] P (x)/Q(x) is summarized in Algorithm 1, and its immediate consequence for computing the N -th term of the linearly recurrent sequence (u n ) n≥0 defined by eq. (1.1) is displayed in Algorithm 2. Algorithm 1 has complexity 2 M(d) log(N + 1) and Algorithm 2 has cost 2 M(d) log(N +1) +M(d), which proves Theorem 1.

Note that Algorithms 1 and 2 use an idea similar to the ones in [START_REF] Bostan | Fast computation of the N th term of an algebraic series over a finite prime field[END_REF][START_REF] Bostan | Fast coefficient computation for algebraic power series in positive characteristic[END_REF] which were dedicated to the larger Algorithm 2 (OneTerm) [START_REF] Householder | Dandelin, Lobačevskiȋ, or Graeffe?[END_REF][START_REF] Pan | Algebraic complexity of computing polynomial zeros[END_REF][START_REF] Pan | Solving a polynomial equation: some history and recent progress[END_REF]. Interestingly, the Graeffe process has been used in a purely algebraic context by Schönhage in [67, §3] for computing the reciprocal of a power series, see also [14, §2]. However, our paper seems to be the first reference where the Graeffe process and the section operators approach are combined together.

Input: rec. (1.1), u 0 , . . . , u d-1 , N Output: u N Assumptions: Γ(x) = x d - d-1 i=0 c i x i with c 0 = 0 1: Q(x) ← x d Γ(1/x) 2: P (x) ← (u 0 + • • • + u d-1 x d-1 ) • Q(x) mod x d 3: return [x N ]P (x)/Q(x)

New algorithm for Fibonacci numbers

To illustrate the mechanism of Algorithm 1, let us instantiate it in the particular case of the Fibonacci sequence

F 0 = 0, F 1 = 1, F n+2 = F n+1 + F n , n ≥ 0. The generating function n≥0 F n x n is x/(1 -x -x 2 ). Therefore, the coefficient F N = [x N ] x 1-x-x 2 is equal to [x N ] x(1 + x -x 2 ) 1 -3x 2 + x 4 = [x N 2 ] x 1-3x+x 2 , if N is even, [x N -1 2 ] 1-x 1-3x+x 2 , else. The computation of F N is reduced to that of a coefficient of the form [x N ] a + bx 1 -cx + x 2 = [x N ] (a + bx)(1 + cx + x 2 ) 1 -(c 2 -2)x 2 + x 4 which is equal to [x N 2 ] a+(bc+a)x 1-(c 2 -2)x+x 2 , if N is even, [x N -1 2 ] (ac+b)+bx 1-(c 2 -2)x+x 2 ,
else. This yields Algorithm 3 for the computation 8 of F N . 8 Notice that the same algorithm can be used to compute efficiently the N -th Fibonacci polynomial, or the N -th Chebyshev polynomial. Fibonacci polynomials in R[t] are defined by

F n+2 (t) = t • F n+1 (t) + Fn(t) with F 0 (t) = 1 and F 1 (t) = 1. It is sufficient to initialize c to t 2 +2 (instead of 3) and b to t when N is even (instead of 0). The complexity of this algorithm is O(M(N )) operations in R, which is quasi-optimal. Algorithm 3 (NewFibo) Input: N Output: F N Assumptions: N ≥ 2 1: c ← 3 2: if N is even then 3: [a, b] ← [0, 1] 4: else 5: [a, b] ← [1, -1] 6: N ← N/2 7: while N > 1 do 8: if N is even then 9: b ← a + b • c 10: else 11: a ← b + a • c 12: c ← c 2 -2 13: N ← N/2 14: return b + a • c
In the direct application of Algorithm 1 to the Fibonacci sequence, the condition N > 1 in line 7 should be N ≥ 1 and b + a • c in line 14 should be a. But, in that case, there is a redundant product c 2 in line 12 when N = 1. This modification saves one multiplication and one subtraction. This saving is crucial over the integers, when we consider bit complexity, since c grows exponentially in the iterations, and the bit complexity of the redundant multiplication c 2 occupies a constant factor of the whole bit complexity.

A close inspection reveals that this algorithm computes F N by a recursive use of the formula

F N = L 2 log N • F N -2 log N + (-1) N • F 2 1+ log N -N ,
which is a particular instance of the classical formula

F n+m = L m F n + (-1) n F m-n
relating the Fibonacci numbers and the Lucas numbers

L n = F n+1 + F n-1 .
When N is a power of 2, then Algorithm 3 degenerates into Algorithm 4. This is equivalent to Algorithm fib(n) in [18, Fig. 6]9 . It uses 2 log(N ) -3 products (of which log(N ) -2 are squarings) and log(N ) -2 subtractions.

When N is arbitrary, Algorithm 3 has essentially the same cost: it uses at most 2 log(N ) -1 products (of which at most log(N )-1 are squarings) and 2 log(N ) -1 additions/subtractions. In contrast, [START_REF] Cull | Computing Fibonacci numbers quickly[END_REF]Fig. 6] uses a more complex algorithm, with higher cost.

A nice feature of Algorithm 3 is not only that it is simple and natural, but also that its arithmetic and bit

Algorithm 4 Input: N

Output: F N Assumptions: N ≥ 2 and N is a power of 2

1: [b, c] ← [1, 3] 2: N ← N/2 3: while N > 2 do 4: b ← b • c 5: c ← c 2 -2 6: N ← N/2 7: return b • c
complexity match the complexities of the state-of-art algorithms for computing Fibonacci numbers [START_REF] Takahashi | A fast algorithm for computing large Fibonacci numbers[END_REF].

The MSB-first algorithm and applications

We present in §3.1 a "most significant bit" (MSB) variant (Algoritm 5) of Algoritm 1. Then we discuss various applications of Algorithms 1 and 5. In §3.2 we design a faster algorithm for polynomial modular exponentiation, that we use in §3.3 to design a faster Fiduccia-like algorithm for computing a slice of d terms of indices N -d+1, . . . , N in ∼ 2 M(d) log N operations.

3.1

The MSB-first algorithm In Fiduccia's algorithm ( §1.3), the N -th coefficient u N in the power series expansion i≥0 u i x i of P/Q is given by the inner product x N mod Γ(x), v 0 , where Γ is the reversal polynomial of Q and v 0 is the vector of initial coefficients

u 0 • • • u d-1 .
Here, x N mod Γ(x) depends only on the linear recurrence equation (1.1), and is independent of the initial terms v 0 . Hence, if we want to compute the N -th terms of k different linearly recurrent sequences that share the same linear recurrence equation (1.1), we can first determine ρ(x) := x N mod Γ(x), and then ρ, v On the other hand, in Algorithm 1 we iteratively update both the denominator and numerator, and each new numerator depends on the original numerator P (x) which encodes the initial d terms of the sequence. Hence, it is not a priori clear how to obtain with Algorithm 1 the good feature of Fiduccia's algorithm.

In this section, we present an algorithm that computes u N with complexity equal to that of Algorithm 1 and which, in addition, achieves the cost O(M(d) log N + kd) for the above problem with k sequences.

While Algorithm 1 looks at N from the least significant bit (LSB), the main algorithm presented in this section (Algoritm 5) looks at N from the most Algorithm 5 (SliceCoeff)

Input: Q(x), N F N,d (1/Q(x)) Assumptions: Q(0) invertible and deg(Q) =: d 1: function SliceCoeff(N , Q(x)) 2:
if N = 0 then 3:

return x d-1 /Q(0) 4: A(x) ← Q(x)Q(-x) A = 2d i=0 A i x i 5: V (x) ← d i=0 A 2i x i 6: W (x) ← SliceCoeff( N/2 , V (x)) 7:
if N is even then 8:

S(x) ← xW (x 2 ) 9:
else 10:

S(x) ← W (x 2 )
11:

B(x) ← Q(-x)S(x) B = 3d-1 i=0 B i x i 12: return d-1 i=0 B d+i x i
significant bit (MSB). In fact, Algoritm 5 is essentially equivalent to "the transposition" of Algorithm 1 in the sense of [START_REF] Bostan | Tellegen's principle into practice[END_REF]. This is the reason why this MSB-first algorithm has the same cost as Algorithm 1. However, in order to keep the presentation self-contained, we are not going to appeal here to the general machinery of algorithmic transposition tools, but rather derive the transposed algorithm "by hand", by a direct reasoning.

In order to compute [x N ] P (x)/Q(x), it is sufficient to compute the (N -d + 1)-th term to the N -th term of 1/Q(x) since the degree of P (x) is at most

d -1. Let F N,d ( i≥0 a i x i ) := d-1 i=0 a N -d+1+i x i . Our goal is to compute F N,d (1/Q(x)
). We have the equalities

F N,d 1 Q(x) = F N,d Q(-x) Q(x)Q(-x) = F N,d Q(-x)x N -2d+1 F N,2d 1 Q(x)Q(-x) = F 2d-1,d Q(-x)F N,2d 1 V (x 2 )
,

where V (x 2 ) := Q(x)Q(-x).
In the second equality, we ignore the terms of 1/V (x 2 ) except for the (N -2d + 1)-th term to the N -th term. In the third equality, we use the fact that

F N,d (xA(x)) = F N -1,d (A(x)). Let now W (x) := F N/2 ,d (1/V (x))
. Then, it is easy to see that

F N,d 1 Q(x) = F 2d-1,d (Q(-x)S(x)) ,
where

S(x) := xW (x 2 ), if N is even W (x 2 ), else.
Algorithm 6 (OneCoeffT)

Input: P (x), Q(x), N Output: [x N ] P (x) Q(x) Assumptions: Q(0) invertible, deg(P ) < deg(Q) =: d 1: U ← F N,d (1/Q(x)) using Algorithm 5 U = u N -d+1 + • • • + u N x d-1 2: return p 0 u N + • • • + p d-1 u N -d+1 P = d-1 i=0 p i x i Algorithm 7 Input: P 1 , . . . , P k , Q, N Output: [x N ] Pj Q , j = 1, . . . , k Assumptions: Q(0) invertible, deg(P ) < deg(Q) =: d 1: U ← F N,d (1/Q(x)) using Algorithm 5 U = u N -d+1 + • • • + u N x d-1 2: return p (j) 0 u N + • • • + p (j) d-1 u N -d+1 , j = 1, . . . , k P j = d-1 i=0 p (j) i x i
The resulting method for computing F N,d (1/Q(x)) is summarized in Algorithm 5, and its applications to the computation of [x N ]P/Q, and to [x N ]P (i) /Q for several i = 1, . . . , k, are displayed in Algorithms 6 and 7.

Let us analyze the complexity of Algorithms 5 and 6 more carefully. At each step, Algorithm 5 computes Q(x)Q(-x) and Q(-x)S(x), where the degrees of Q(x) and S(x) are d and at most 2d -1, respectively. Hence a direct analysis concludes that its complexity is 3 M(d) log N operations in R. However, an improvement comes from the remark that not all coefficients of Q(-x)S(x) are needed: it is sufficient to compute the dth coefficient to the (2d-1)-th coefficient of Q(-x)S(x). This operation is known as "the middle product", and can be performed with the same arithmetic complexity as the standard product of two polynomials of degree d, plus a few d additional operations [START_REF] Hanrot | The middle product algorithm[END_REF][START_REF] Bostan | Tellegen's principle into practice[END_REF]. Therefore, if steps 11 and 12 of Algorithm 5 are performed "at once" using a middle product, then the arithmetic complexity drops to (2 M(d) + d) log N . This complexity is also inherited by Algorithm 6, which uses at most 2d additional operations in the last step.

It should be obvious at this point that the slight variant Algorithm 7 of Algorithm 6 achieves arithmetic complexity O(M(d) log N + kd) for the aforementioned problem with k sequences, and more precisely its cost is of at most (2 M(d) + d) log N + 2kd operations in R.

In conclusion, Algorithm 5 achieves the same arithmetic complexity as Algorithm 1 and it extends to Algorithms 6 and 7. All algorithmic techniques specific to the FFT setting, that we will describe in Section 4, can also be applied to Algorithms 5, 6 and 7, yielding the same complexity gains.

(NewModExp) Input: Γ(x), N

Output: x N mod Γ(x) Assumptions: lcoeff(Γ) invertible, Γ(0) = 0 and deg(Γ) =: d

1: Q(x) ← x d Γ(1/x) 2: u(x) ← F N,d (1/Q(x)) using Algorithm 5 3: v(x) ← u(x)Q(x) mod x d 4: return v(1/x)x d-1

Faster modular exponentiation

The algorithms of §3.1 are not only well-suited to compute the N -th terms of several sequences satisfying the same recurrence. In this section, we show that they also permit a surprising application to the computation of polynomial modular exponentiations. This fact has many consequences, since modular exponentiation is a central algorithmic task in algebraic computations. In §3.3, we will discuss a first application in relation with the main topic of our article. Namely, we will design a new Fiduccia-style algorithm for the computation of the Nth term, and actually of a whole slice of

k ≥ d terms, in (2 M(d) + d) log N + O((k + d)M(d)/d) ring operations.
More consequences will be separately discussed in §3. [START_REF] Bernstein | Removing redundancy in highprecision Newton iteration[END_REF].

Assume we are given a polynomial Γ(x) ∈ R[x] of degree d, an integer N , and that we want to compute ρ(x) := x N mod Γ(x). Without loss of generality, we may assume Γ(0) = 0. Let Q(x) ∈ R[x] be the reversal of Γ(x), that is Q(x) := x d Γ(1/x). Let us denote the power series expansion of 1/Q by i≥0 a i x i . Then, equation (1.3) implies that

(3.5) a N • • • a N +d-1 = r × H, where r = r 0 • • • r d-1 with ρ = d-1 i=0 r i x i and H is the Hankel matrix H :=      a 0 • • • a d-1 a 1 • • • a d . . . a d-1 • • • a 2d-2      . Note that the matrix H is invertible, as its determinant is equal (up to a sign) to ([x d ]Q) d-1 = Γ(0) d-1 .
Therefore, r (and thus ρ) can be found by (1) computing a N • • • a N +d-1 using Algorithm 5;

(2) solving the Hankel linear system (3.5).

Step (1) has arithmetic cost (2 M(d) + d) log(N + d), while step (2) has negligible cost O(M(d) log d) using [START_REF] Brent | Fast solution of Toeplitz systems of equations and computation of Padé approximants[END_REF], see also [START_REF] Bini | Polynomial and matrix computations[END_REF]Ch. 2,[START_REF] Bernstein | Fast multiplication and its applications[END_REF].

It is actually possible to improve a bit more on this algorithm, by using the next lemma.

Lemma 2. Let N ∈ N and let

Γ(x) ∈ R[x] be of degree d with Γ(0) = 0. Let Q(x) ∈ R[x] be its reversal, Q(x) := x d Γ(1/x). Denote its reciprocal 1/Q by i≥0 a i x i , and let u(x) be F N,d (1/Q(x)) = a N -d+1 + • • •+a N x d-1 . Define v(x) to be u(x)Q(x) mod x d . Then x N mod Γ(x) = v(1/x)x d-1 . Proof. Write the Euclidean division x N = L(x) • Γ(x) + ρ(x), where deg(L) = N -d and ρ(x) = r 0 + • • • + r d-1 x d-1
. Replacing x by 1/x on both sides, and then multiplying by

x N yields 1 = L rev (x) • Q(x) + x N -d+1 • ρ(x), where L rev (x) = x N -d • L(1/x) and ρ(x) = x d-1 • ρ(1/x). In other words 1 Q(x) = L rev (x) + x N -d+1 • ρ(x) Q(x) .
Since L rev (x) has degree at most N -d, it follows that u(x) = ρ(x) Q(x) mod x d . Therefore, ρ(x) is equal to v(x), and the conclusion follows.

The merit of Lemma 2 is that it shows that computing x N mod Γ(x) can be reduced to computing F N,d (1/Q(x)), plus a few additional operations with negligible cost M(d).

The resulting method is presented in Algorithm 8, whose complexity is

(2 M(d) + d) log N + M(d). This proves Theorem 4.
Note that Algorithm 8 is simpler, and faster by a factor of 1.5, than the classical algorithm based on binary powering in the quotient ring R[x]/(Γ(x)). Algorithm 5, and hence also Algorithm 8, admits a specialization into the FFT setting (Algorithm 13 in §4) with complexity ∼ 2 3 M(d) log N . Similarly to the case of Algorithm 11 in §4, this FFT variant of Algorithm 8 is faster by a factor of 2.5 than Shoup's (comparatively simple) algorithm [70, §7.3], and by a factor of 1.625 than the (much more complex) algorithm of Mihăilescu [START_REF] Mihăilescu | Fast convolutions meet Montgomery[END_REF].

This speed-up might be beneficial for instance in applications to polynomial factoring in F p [x], where one time-consuming step to factor f ∈ F p [x] is the computation of x p mod f , see [START_REF] Zur Gathen | Modern computer algebra[END_REF]Algorithms 14.3,14.8,14.13,14.15,14.31,14.33 and 14.36], and also [START_REF] Shoup | A new polynomial factorization algorithm and its implementation[END_REF][START_REF] Lecerf | New recombination algorithms for bivariate polynomial factorization based on Hensel lifting[END_REF].

It is also so in point-counting methods such as Schoof's algorithm and the Schoof-Elkies-Atkin (SEA) algorithm [7, Ch. VII], the second one being the best known method for counting the number of points of elliptic curves defined over finite fields of large characteristic. Indeed, the bulks of these algorithms are computations of x q modulo the "division polynomial" f (x) and of x q modulo the "modular polynomial" Φ (x), where = O(log(q)) and deg(f

) = O( 2 ), deg(Φ ) = O( ).
As a final remark, note that while Fiduccia's algorithm shows that computing the terms of indices

Assumptions: Γ(x) = x d - d-1 i=0 c i x i with c 0 = 0 1: ρ(x) ← x N mod Γ(x) using Algorithm 8 2: U (x) ← u 0 + • • • + u 2d-2 x 2d-2
using Algorithm in [69, p. 18]

3: V (x) ← U (x) • (x d • ρ(1/x)) V = d-1 i=0 v i x i 4: return [v d , . . . , v 2d-1 ]
N, . . . , N + d -1 of a linearly recurrent sequence of order d can be reduced to polynomial modular exponentiation (x N mod Γ(x)), Algorithm 8 shows that the converse is also true: polynomial modular exponentiation can be reduced to computing the terms of indices N, . . . , N + d -1 of a linearly recurrent sequence of order d. Therefore, these two problems are computationally equivalent. To our knowledge, this important fact seems not to have been noticed before.

A new Fiduccia-style algorithm

We conclude this section by discussing a straightforward application of Algorithm 8. This is based on the next equality, generalizing (3.5) to any k ≥ 1:

(3.6) u N • • • u N +k-1 = r × H k ,
where as before

r = r 0 • • • r d-1 is the coefficients vector of ρ = d-1 i=0 r i x i , with ρ = x N mod Γ(x)
, and H k is the Hankel matrix

H k :=      u 0 • • • u d-1 • • • • • • u k-1 u 1 • • • u d • • • • • • u k . . . . . . . . . . . . u d-1 • • • u 2d-2 • • • • • • • • • u k+d-2      .
The matrix H k is built upon the first terms of the sequence (u n ) n≥0 satisfying recurrence (1.1) with characteristic polynomial Γ = x dd-1 i=0 c i x i , or equivalently, from the power series expansion of the rational function

P/Q with Q(x) = x d Γ(1/x).
Note that the entries of H k can be computed either from P and Q, or from the recurrence (1.1) together with the initial terms u 0 , . . . , u d-1 , using O((k + d) M(d)/d) arithmetic operations, by the algorithm in [START_REF] Shoup | A fast deterministic algorithm for factoring polynomials over finite fields of small characteristic[END_REF]Thm. 3.1], see also [8, §5]. To compute u N , . . . , u N +k-1 it thus only remains to perform the vector-matrix product (3.6).

When k = 1, the product r×H 1 costs 2d operations and it yields the term u N .

When k ≥ d, the product r × H k can be reduced to the polynomial multiplication of

r d-1 + • • • + r 0 x d-1 by k+d-2 i=0
u i x i , and this can be performed using k+d d M(d) arithmetic operations. As a consequence, the whole slice of coefficients u N +i = [x N +i ]P/Q for i = 0, . . . , k -1, can be computed using Algorithm 8 and eq. (3.6) for a total cost of arithmetic operations of

(2 M(d) + d) log N + O k + d d M(d) .
When k = d, this proves Theorem 3. The corresponding algorithm is given as Algorithm 9.

We emphasize that this variant of Fiduccia's algorithm is different from Algorithm 1. It is actually a bit slower than Algorithm 1 when k = 1. However, when k > 1 terms are to be computed, it should be preferred to repeating k times Algorithm 1. It also compares favorably with Fiduccia's original algorithm, whose adaption to k terms has arithmetic complexity

3 M(d) log N + O d log(N ) + k + d d M(d) .

Applications

In this section, we discuss more applications of the MSB-first algorithms (Algorithm 5 and 8) presented in §3.1 and §3.2. We deal with the case of multiplicities ( §3.4.1), and explain a new way to speed up computations in that case. Then, we address another application, to faster powering of matrices ( §3.4.2).

To simplify matters, we assume in this section that R = K is a field.

3.4.1

The case with multiplicities Hyun and his co-authors [START_REF] Hyun | A fast algorithm for solving linearly recurrent sequences[END_REF][START_REF] Hyun | Change of basis for m-primary ideals in one and two variables[END_REF] addressed the following question: is it possible to compute faster the N -th term of a linearly recurrent sequence when the characteristic polynomial of the recurrence has multiple roots? By the Chinese Remainder Theorem, it is sufficient to focus on the case where the characteristic polynomial is a pure power of a squarefree polynomial. In other words, the main step of [42, Algorithm 1] is to compute x N mod Q, where Q = (Q ) m and Q is the squarefree part of Q. Under suitable invertibility conditions, the problem is solved in [START_REF] Hyun | A fast algorithm for solving linearly recurrent sequences[END_REF][START_REF] Hyun | Change of basis for m-primary ideals in one and two variables[END_REF] in

O(M(d ) log N + M(d) log d) operations in K, where d = deg(Q ) and d = deg(Q) = m • d .
This cost is obtained using an algorithm based on bivariate computations, using the isomorphisms between K[x]/(Q) and K[y, x]/(Q (y), (x -y) m ) made effective by the so-called tangling / untangling operations. We now propose an alternatively fast, but simpler, algorithm with the same cost.

Let us explain this on an example, for "multiple-Fibonacci numbers", that is when

Q has the form (Q ) m , with Q = 1 -x -x 2 and d = 2, d = 2m.
Assume we want to compute the N -th coefficient u N in the power series expansion of (x/(1 -x -x 2 )) m . cost of Fiduccia's algorithm, and also of our new algorithms, is O(M(m) • log N ).

Let us explain how we can lower this to O(log N + M(m) log m). The starting point is the observation that, by the structure theorem of linearly recurrent sequences [15, §A.(I)] (see also [62, §2]), u N is of the form w m (N )φ N +v m (N )ψ N , where φ and ψ are the two roots of 1+x = x 2 and w m , v m are polynomials in K[x] of degree less than m. By an easy liner algebra argument,

u N is thus equal to W m (N )F N + V m (N )F N +1
, where W m (x) and V m (x) are polynomials in K[x] of degree less than m. These polynomials can be computed by (structured) linear algebra from the first 2m values of the sequence (u n ), in complexity O(M(m) log m).

For instance, when d = 2, we have U 2 (x) = -(x + 1)/5 and V 2 (x) = 2x/5. Once U m and V m are determined, it remains to compute F N and F N +1 using Algorithm 9 in O(log N ) operations in K, then to return the value W m (N )

• F N + V m (N ) • F N +1 .
The arguments extend to the general case and yields an algorithm of arithmetic complexity

(2 M(d ) + d ) log N + O(M(d) log d).

Faster powering of matrices

Assume we are given a matrix M ∈ M d (K), an integer N , and that we want to compute the N -th power M N of M .

The arithmetic complexity of binary powering in M d (K) is O(d θ log N ) operations in K, where as before θ ∈ [2, 3] is any feasible exponent for matrix multiplication in M d (K).

A better algorithm consists in first computing the characteristic polynomial Γ(x) of the matrix M , then the remainder ρ(x) := x N mod Γ(x), and finally evaluating the polynomial ρ(x) at M . By the Cayley-Hamilton theorem, ρ(M ) = M N . The most costly step is the computation of ρ, which can be done as explained in §3.2 using ∼ 2 M(d) log(N ) operations in K. The cost of the other two steps is independent of N , and it is respectively O(d θ log d) [START_REF] Keller-Gehrig | Fast algorithms for the characteristic polynomial[END_REF] and O(d θ+ 1 2 ), this last cost being achieved using the Paterson-Stockmeyer baby-step/giant-step algorithm [START_REF] Paterson | On the number of nonscalar multiplications necessary to evaluate polynomials[END_REF]. The total cost of this algorithm is (2

M(d) + d) log(N ) + O(d θ+ 1
2 ). Note that a faster variant (w.r.t. d), of cost (2 M(d) + d) log(N ) + O(d θ log d), can be obtained using [START_REF] Giesbrecht | Nearly optimal algorithms for canonical matrix forms[END_REF]Corollary 7.4]. The corresponding algorithm is based on the computation of the Frobenius (blockcompanion) form of the matrix M , followed by the powering of companion matrices, which again reduces to modular exponentiation.

Analysis under the FFT multiplication model

In this section, we specialize, optimize and analyze the generic Algorithm 1 to the FFT setting, in which polynomial products are assumed to be performed using the discrete Fourier transform (DFT), and its inverse.

To do this, we will assume that the base ring R possesses roots of unity of sufficiently high order. To simplify the exposition, R will be supposed to be a field, but the arguments also apply without this assumption, modulo some technical complications, see [27, §8.2].

Discrete Fourier Transform for polynomial products

Let K be a field with a primitive n-th root ω n of unity. Let A ∈ K[x] be a polynomial of degree at most d ≤ n -1. The DFT A of A is defined by

A y := A(ω -y n ) = n-1 i=0 A i ω -yi n for y = 0, 1, . . . , n -1.
Here, A y = 0 for y > d. It is classical that the DFT map is an invertible K-linear transform from K n to itself, and that the polynomial A can be retrieved from its DFT A using the formulas

A i = 1 n n-1 y=0 A y ω yi n for i = 0, 1, . . . , n -1.
For computing the polynomial multiplication

C(x) = A(x)B(x) for given A(x), B(x) ∈ K[x] of degree at most d, it is sufficient to compute the DFT of C(x) for n ≥ 2d + 1. Since C y = C(ω -y n ) = A(ω -y n )B(ω -y n ) = A y B y
, the polynomial C(x) can be computed using two DFTs and one inverse DFT.

Let E(n) be an arithmetic complexity for computing a DFT of length n. Then the cost of polynomial multiplication in K[x] is governed by

M(d) = 3 E(2d) + O(d).

Fast Fourier Transform

In this subsection, we briefly recall the Fast Fourier Transform (FFT), which gives the quasi-linear estimate

E(n) = O(n log n).
Assume n is even. Then, for y = 0, 1, . . . , n/2 -1 we have Thus, computing a DFT in size n reduces to two DFTs in size n/2. More precisely, E(n) ≤ 2 E(n/2)+(3/2)n. If n is a power of two, n = 2 k , and if the field K contains a primitive 2 k -th root of unity (as is the case for instance when K = C, K = F p for a prime number p satisfying 2 k | p -1), this reduction can be repeated k = log n times, and it yields the estimate E(n) = 3 2 n log n. The corresponding algorithm is called the decimationin-time Cooley-Tukey fast Fourier transform [START_REF] Cooley | An algorithm for the machine calculation of complex Fourier series[END_REF], see also [5, §2].

A y = n/2-1 i=0 A 2i ω -y(2i) n + n/2-1 i=0 A 2i+1 ω -y(2i+1) n = n/2-1 i=0 A 2i ω -yi n/2 + ω -y n • n/2-1 i=0 A 2i+1 ω -yi
By the arguments of §4.1, we conclude that polynomial multiplication in K[x] can be performed in arithmetic complexity

M(d) = 9 d log d + O(d).

Efficiently doubling the length of a DFT

In the FFT setting, it is useful for many applications to compute efficiently a DFT of length 2n starting from DFT of length n.

Assume n ≥ d + 1 and we have at our disposal the DFT A of A, of length n. Assume that we want to compute the DFT A (2n) of length 2n.

The simplest algorithm is to apply the inverse DFT of length n to obtain A, and then to apply the DFT of length 2n to A. This costs E(n) + E(2n) arithmetic operations, that is 9 2 n log n + 3n operations in K. This algorithm can be improved using the following formulas

A (2n) 2y = 2n-1 i=0 A i ω -2yi 2n = n-1 i=0 A i ω -yi n = A y , A (2n) 2y+1 = 2n-1 i=0 A i ω -(2y+1)i 2n = n-1 i=0 ω -i 2n A i ω -yi n = B y ,
where B i := ω -i 2n A i for i = 0, 1, . . . , n -1. We obtain Algorithm 10 with arithmetic complexity 2 E(n)+n, i.e. 3n log n+n, [5, §12.8], see also [START_REF] Bernstein | Removing redundancy in highprecision Newton iteration[END_REF][START_REF] Mezzarobba | NumGfun: a package for numerical and analytic computation and D-finite functions[END_REF] 10 . Compared with the direct algorithm, the gain is roughly a factor of 3/2.

Algorithm 1 in the FFT setting

Recall that our main objective is, given P, Q in K[x] with d = 10 This "FFT doubling" trick is sometimes attributed to R. Kramer ( 2004), but we were not able to locate Kramer's paper. Let k be the minimum integer satisfying 2 k ≥ 2d+1. Assume that there exists a primitive 2 k -th root of unity in K. In this case, we can employ an FFT-based polynomial multiplication in K[x]. In each iteration of Algorithm 1, it is sufficient to compute P (x)Q(-x) and Q(x)Q(-x). Here, only two FFTs and two inverse FFT of length 2 k are needed since

Algorithm 10 Input: DFT

n (A) Output: DFT 2n (A) 1: function DoubleDFT( A) 2: A ← IDFT n ( A) 3: B i ← ω -i 2n A i for i = 0, 1, . . . , n -1 4: B ← DFT n (B) 5
Q - y = Q ȳ for Q -(x) := Q(-x) where ȳ := y + 2 k-1 if y < 2 k-1 and ȳ := y -2 k-1 if y ≥ 2 k-1 . Hence, the arithmetic complexity S(d) for a single step in Algorithm 1 satisfies S(d) ≤ 4 E(2 k ) + O(2 k ).
In the following we will show the improved estimate

S(d) ≤ 4 E(2 k-1 ) + O(2 k ).
Before entering the while loop in Algorithm 1, the DFTs P and Q of P (x) and Q(x) of length 2 k are computed, respectively. Inside the while loop, P and Q are updated. The recursive formula (4.7) for the decimation-in-time Cooley-Tukey FFT is equivalent to

A e y A o y = 1 2 1 0 0 ω y n 1 1 1 -1 A y A 2 k-1 +y for y = 0, 1, . . . , 2 k-1 .
By using this formula, A e (or A o ) can be computed with O(2 k ) operations from A. By using Algorithm 10, we obtain the updated P from A e or A o . The algorithm is summarized in Algorithm 11. In each step, DoubleDFT is called twice. Hence, the total arithmetic complexity of Algorithm 11 is

(4 E(2 k-1 ) + O(2 k )) • log N.
When d is of the form 2 -1, then11 one can take k = + 1 and the cost simplifies to

T(N, d) = 4 E(d) log N + O(d log N ), 11 (OneCoeff-FFT) Input: P (x), Q(x), N Output: [x N ] P (x) Q(x) 1: P ← DFT 2 k (P ) 2: Q ← DFT 2 k (Q) 3: while N ≥ 1 do 4: U y ← P y Q ȳ for y = 0, 1, . . . , 2 k -1 5: if N is even then 6: U e y ← ( U y + U y+2 k-1 )/2 for y = 0, 1, . . . , 2 k-1 -1 7: P ← DoubleDFT( U e y ) 8: else 9: U o y ← ω y 2 k ( U y -U y+2 k-1 )/2 for y = 0, 1, . . . , 2 k-1 -1 10: P ← DoubleDFT( U o y ) 11: A y ← Q y Q ȳ for y = 0, 1, . . . , 2 k-1 -1 12: Q ← DoubleDFT( A) 13: N ← N/2 14: P (0) ← 2 k -1 y=0 P y 15: Q(0) ← 2 k -1 y=0 Q y 16: return P (0)/Q(0) or, equivalently T(N, d) = 6 d log d log N + O(d log N ).
The (striking) conclusion of this analysis is that, in the FFT setting, our (variant of the) algorithm for computing the N -th term of P/Q uses much less operations than in the general case, namely Note that the complexity bound (4.8) compares favorably with Fidducia's algorithm combined with the best algorithms for modular squaring. For instance, Shoup's algorithm [70, §7.3] computes one modular squaring in the FFT setting using ∼ 5 3 M(d) arithmetic operations, while Mihăilescu's algorithm [53, Table 1] (based on Montgomery's algorithm [START_REF] Montgomery | Modular multiplication without trial division[END_REF]) uses roughly ∼ 13 12 M(d) arithmetic operations. Our bound (4.8) is better by a factor of 2.5 than Shoup's (comparatively simple) algorithm, and by a factor of 1.625 than the (much more complex) algorithm by Mihăilescu. Let us point out that all the other algorithms admit similarly fast versions in the FFT setting. We will however not give them in full detail here, mainly for space reasons. x) and one middle-product Q(-x)S(x). The standard product can be performed by two FFTs of length 2 k-1 as in Algorithm 11. For the middle-product, we first compute the cyclic product. i.e., the product in K[x]/(x 2 k -1), between Q(-x) and S(x), and then extract the second half of the cyclic product. For computing the DFT of the second half A sec (x) :=

2 k (A) Output: DFT 2 k-1 (A sec ) 1: function SecondHalfDFT( A) 2: B y ← A 2y+1 for y = 0, 1, . . . , 2 k-1 -1 3: B ← IDFT 2 k-1 ( B) 4: C i ← ω i 2 k B i for i = 0, 1, . . . , 2 k-1 -1 5: C ← DFT 2 k-1 (C) 6: D y ← ( A 2y -C y )/2 for y = 0, 1, . . . , 2 k-1 -1 7: return D arithmetic complexity ((2/3)M(d) + O(d)) log N . Algo- rithm 5 includes one standard product Q(x)Q(-
2 k-1 -1 i=0
A i+2 k-1 x i , we use the formula for the decimation-in-frequency FFT.

A 2y = 2 k -1 i=0 ω -2iy 2 k A i = 2 k-1 -1 i=0 ω -iy 2 k-1 (A i + A i+2 k-1 ) A 2y+1 = 2 k -1 i=0 ω -i(2y+1) 2 k A i = 2 k-1 -1 i=0 ω -iy 2 k-1 ω -i 2 k (A i -A i+2 k-1 ).
This formula gives Algorithm 12 which transforms A to A sec . The arithmetic complexity of Algorithm 12 is 2E(2 k-1 ) + O(2 k ). The whole algorithm is shown in Algorithm 13. The input of the algorithm is N and Q that is the DFT of Q of length 2 k . The output of the algorithm is the DFT of F N,2 k-1 (1/Q(x)) of length 2 k-1 rather than length 2 k since the DFT of W (x 2 ) of length 2 k can be obtained efficiently from the DFT of W (x) of length 2 k-1 . Algorithm 13 calls one DoubleDFT and one SecondHalfDFT. Hence, similarly to Algorithm 11, the arithmetic complexity of Algorithm 13 is (4E(2 k-1 ) + O(2 k )) log N .

Conclusion

We have proposed several algorithmic contributions to the classical field of linearly recurrent sequences. Firstly, we have designed a simple and fast algorithm for computing the N -th term of a linearly recurrent sequence of order d, using ∼ 2 M(d) log N arithmetic operations, which is faster by a factor of 1.5 than the state-of-the-art 1985 algorithm due to Fiduccia [START_REF] Fiduccia | An efficient formula for linear recurrences[END_REF]. When combined with FFT techniques, the algorithm has even better arithmetic complexity ∼ if N = 0 then 3:

Q(0) ← 2 k -1 y=0 Q y 4: return [1 ω 2 k ω 2 2 k • • • ω 2 k -1 2 k
]/Q(0)

5:

A y ← Q y Q ȳ for y = 0, 1, . . . , 2 k-1 -1 6:

V ← DoubleDFT( A)

7:

W ← SliceCoeff-FFT( N/2 , V )

8:
if N is even then B y ← S y Q ȳ for y = 0, 1, . . . , 2 k -1

13: return SecondHalfDFT( B)
which is faster than the fastest variant of Fiduccia's algorithm in the FFT setting by a factor of 1.625. The new algorithms are based on a new method (Algorithm 1) for computing the N -th coefficient of a rational power series. Secondly, using algorithmic transposition techniques, we have derived from Algorithm 1 a new method (Algorithm 5) for computing simultaneously the coefficients of indices N -d + 1, . . . , N in the power series expansion of the reciprocal of a degree-d polynomial, using again ∼ 2 M(d) log N arithmetic operations. Using Algorithm 5, we have designed a new algorithm for computing the remainder of x N modulo a given polynomial of degree d, using ∼ 2 M(d) log N arithmetic operations as well. This is better by a factor of 1.5 than the previous best algorithm for modular exponentiation, with an even better speed-up in the FFT setting, as for Algorithm 1. Combined with the basic idea of Fiduccia's algorithm, our new algorithm for modular exponentiation yields a faster Fiduccia-like algorithm (by the aforementioned constant factors) that computes a slice of d consecutive terms (of indices N -d + 1, . . . , N ) of a linearly recurrent sequence of order d using ∼ 2 M(d) log N arithmetic operations.

Thirdly, we have discussed applications of the new algorithms to a few other algorithmic problems, including powering of matrices and the computation of terms of linearly recurrent sequences when the recurrence has roots with (high) multiplicities.

As future work, we plan to investigate further the full power of our technique. To which extent can it be generalized to larger classes of power series? For instance, although it perfectly works for bivariate rational power series U (x, y), the corresponding method does not directly provide a O(log N )-algorithm for computing the (N, N )-th coefficient u N,N , the reason being that the log N new bivariate recurrences produced by the Graeffe process do not have constant orders, as in the univariate case. This is disappointing, but after all not surprising, because the generating function of the sequence (u n,n ) n is known to be algebraic, but not rational anymore [START_REF] Pólya | Sur les séries entières, dont la somme est une fonction algébrique[END_REF]. As of today, no algorithm is known for computing the N -th coefficient of an algebraic power series faster than in the P-recursive case (P), namely in a number of ring operations almost linear in √ N .

  at the cost of O(log N ) multiplications in A. Each multiplication in A may be performed using O(M(d)) operations in R [27, Ch. 9, Corollary 9.7], where M(d) denotes the arithmetic cost of polynomial multiplication in R[x] in degree d. Using Fast Fourier Transform (FFT) methods, one may take M(d) = O(d log d) when R contains enough roots of unity, and M(d) = O(d log d log log d) in general [27, Ch. 8].

  i = 1, . . . , k, where v (i) 0 denotes the vector of d initial terms of the i-th sequence. The total arithmetic complexity of this algorithm is O(M(d) log N + kd); this is faster than Fiduccia's algorithm repeated independently k times, with cost O(k M(d) log N ).

Algorithm 9

 9 Input: rec. (1.1), u , . . . , u d-1 , N Output: u N , . . . , u N +d-1

n/ 2 =

 2 A e y + ω -y n A o y for A e (x) := n/2 i=0 A 2i x i and A o (x) := n/2 i=0 A 2i+1 x i . Similarly, we have A n/2+y = A e y -ω -y n A o y . We therefore obtain the following matrix equation y A y = 0, 1, . . . , n/2 -1.

7 :

 7 y for y = 0, 1, . . . , n -1 y for y = 0, 1, . . . , n -1 return A(2n) deg(Q) > deg(P ), to compute the N -th coefficient u N in the series expansion of P/Q.

  log N, while for a generic multiplication algorithm the cost is ∼ 2 M(d) log N . This proves Theorem 2.

4. 5

 5 Algorithm 5 in the FFT setting In this section, we present Algorithm 5 in the FFT setting with Algorithm 12 Input: DFT

2 3 M 1 :

 31 (d) log N Algorithm 13 (SliceCoeff-FFT) Input: DFT 2 (Q), N Output: DFT 2 k-1 (F N,2 k-1 (1/Q(x))) function SliceCoeff-FFT(N , Q) 2:

9 : 1 10: else 11 :

 9111 S y ← ω -y 2 k W y mod 2 k-1 for y = 0, 1, . . . , 2 k -S y ← W y mod 2 k-1 for y = 0, 1, . . . , 2 k -1 12:

The idea already appears in the 1982 conference paper[START_REF] Fiduccia | The n-th power of a companion matrix: Fast solutions to linear recurrences[END_REF]. We have discovered that the same algorithm had been sketched by D. Knuth in the corrections/changes to[START_REF] Knuth | The art of computer programming[END_REF] published in 1981 in[47, p. 28], where he attributes the result to R. Brent. Almost surely, C. Fiduccia was not aware about this fact.

Note that multiplying by x in A is much easier and has linear arithmetic cost O(d).

Here, and in what follows, the "equivalent" sign ∼ is to be understood for fixed N and d going to infinity.

This algorithm had also appeared before, in Knuth's book [45, p. 552, second solution].

In the general case, it might be useful to use the Truncated Fourier Transform (TFT), which smoothes the "jumps" in complexity exhibited by FFT algorithms[START_REF] Van Der Hoeven | The truncated Fourier transform and applications[END_REF][START_REF] Harvey | An in-place truncated Fourier transform and applications to polynomial multiplication[END_REF][START_REF] Arnold | A new truncated Fourier transform algorithm[END_REF].

Acknowledgements. We are grateful to the three reviewers for their kind and constructive remarks. Our special thanks go to Kevin Atienza, whose editorial on https://discuss.codechef.com was our initial source of inspiration, and to Bruno Salvy and Sergey Yurkevich, for their careful reading of a first draft of this work. A. Bostan was supported in part by DeRerumNatura ANR-19-CE40-0018. R. Mori was supported in part by JST PRESTO Grant #JPMJPR1867 and JSPS KAKENHI Grant #JP17K17711, #JP18H04090 and #JP20H04138.