Alin Bostan
email: alin.bostan@inria.fr

Sergey Yurkevich
email: sergey.yurkevich@univie.ac.at

Fast Computation of the N -th Term of a q-Holonomic Sequence and Applications

In 1977, Strassen invented a famous baby-step/giant-step algorithm that computes the factorial N ! in arithmetic complexity quasi-linear in √ N . In 1988, the Chudnovsky brothers generalized Strassen's algorithm to the computation of the N -th term of any holonomic sequence in essentially the same arithmetic complexity. We design q-analogues of these algorithms. We first extend Strassen's algorithm to the computation of the qfactorial of N , then Chudnovskys' algorithm to the computation of the N -th term of any q-holonomic sequence. Both algorithms work in arithmetic complexity quasi-linear in √ N ; surprisingly, they are simpler than their analogues in the holonomic case. We provide a detailed cost analysis, in both arithmetic and bit complexity models. Moreover, we describe various algorithmic consequences, including the acceleration of polynomial and rational solving of linear q-differential equations, and the fast evaluation of large classes of polynomials, including a family recently considered by Nogneng and Schost.

Introduction

A classical question in algebraic complexity theory is: how fast can one evaluate a univariate polynomial at one point? The precise formulation of this question depends on the model of computation. We will mainly focus on the arithmetic complexity model, in which one counts base field operations at unit cost.

Horner's rule evaluates a polynomial P in O(deg(P )) operations. Ostrowski [START_REF] Ostrowski | On two problems in abstract algebra connected with Horner's rule[END_REF] conjectured in 1954 that this is optimal for generic polynomials, i.e., whose coefficients are algebraically independent over the prime subfield. This optimality result was proved a few years later by Pan [START_REF] Pan | Methods of computing values of polynomials[END_REF].

However, most polynomials that one might wish to evaluate "in practice" have coefficients which are not algebraically independent. Paterson and Stockmeyer [START_REF] Paterson | On the number of nonscalar multiplications necessary to evaluate polynomials[END_REF] showed, using the baby-step/giant-step technique, that for any field K, an arbitrary polynomial P ∈ K[x] of degree N can be evaluated at any point in an arbitrary K-algebra A using O( √ N ) nonscalar multiplications, i.e., multiplications in A. However, their algorithm uses a linear amount of scalar multiplications, so it is not well adapted to the evaluation at points from the base field K, since in this case the total arithmetic complexity, counted in terms of operations in K, remains linear in N .

For some families of polynomials, one can do much better. Typical examples are x N and P N (x) :

= x N -1 + • • • + x + 1,
which can be evaluated by the square-and-multiply technique in O(log N ) operations.

(Note that for P N (x) such a fast algorithm needs to perform division.) By contrast, a family (F n (x)) n of univariate polynomials is called hard to compute if for large enough N , the complexity of the evaluation of F N grows at least like a power in deg(F N ), whatever the algorithm used.

Paterson and Stockmeyer [START_REF] Paterson | Bounds on the evaluation time for rational polynomial[END_REF][START_REF] Paterson | On the number of nonscalar multiplications necessary to evaluate polynomials[END_REF] proved the existence of polynomials in Q[x] which are hard to compute (note that this does not follow from Pan's result [START_REF] Pan | Methods of computing values of polynomials[END_REF]). However, their proof was based on a non-constructive argument. Specific families of hard-tocompute polynomials were first exhibited by Strassen [START_REF] Strassen | Polynomials with rational coefficients which are hard to compute[END_REF]. For instance, he proved that for large N , the polynomial N ℓ=0 2 2 ℓ x ℓ needs at least N/(3 log N ) operations to be evaluated. The techniques were refined and improved by Borodin and Cook [START_REF] Borodin | On the number of additions to compute specific polynomials[END_REF], Lipton [START_REF] Lipton | Polynomials with 0 -1 coefficients that are hard to evaluate[END_REF] and Schnorr [START_REF] Schnorr | Improved lower bounds on the number of multiplications / divisions which are necessary to evaluate polynomials[END_REF], who produced explicit examples of degree-N polynomials whose evaluation requires a number of operations linear in √ N . Subsequently, various methods have been developed to produce similar results on lower bounds, e.g., by Heintz and Sieveking [START_REF] Heintz | Lower bounds for polynomials with algebraic coefficients[END_REF] using algebraic geometry, and by Aldaz et al. [START_REF] Aldaz | A new method to obtain lower bounds for polynomial evaluation[END_REF] using a combinatorial approach. The topic is vast and very well summarized in the book by Bürgisser, Clausen and Shokrollahi [START_REF] Bürgisser | Algebraic complexity theory[END_REF].

In this article, we focus on upper bounds, that is on the design of fast algorithms for special families of polynomials, which are hard to compute, but easier to evaluate than generic polynomials. For instance, for the degree- N 2 polynomial

Q N (x) := P 1 (x) • • • P N (x),
a complexity in O(N ) is clearly achievable. We will see in §2.1 that one can do better, and attain a cost which is almost linear in √ N (up to logarithmic factors in N ). Another striking example is

R N (x) := N ℓ=0
x ℓ 2 , of degree N 2 , and whose evaluation can also be performed in complexity quasi-linear in √ N , as shown recently by Nogneng and Schost [START_REF] Nogneng | On the evaluation of some sparse polynomials[END_REF] (see §2.2). In both cases, these complexities are obtained by clever although somehow ad-hoc algorithms. The starting point of our work was the question whether these algorithms for Q N (x) and R N (x) could be treated in a unified way, which would allow to evaluate other families of polynomials in a similar complexity.

The answer to this question turns out to be positive. The key idea, very simple and natural, is to view both examples as particular cases of the following general question: Given a q-holonomic sequence, that is, a sequence satisfying a linear recurrence with polynomial coefficients in q and q n , how fast can one compute its N -th term?

In the more classical case of holonomic sequences (satisfying linear recurrences with polynomial coefficients in the index n), fast algorithms exist for the computation of the N -th term. They rely on a basic block, which is the computation of the factorial term N ! in arithmetic complexity quasi-linear in √ N , using an algorithm due to Strassen [START_REF] Strassen | Einige Resultate über Berechnungskomplexität[END_REF]. The Chudnovsky brothers extended in [START_REF] Chudnovsky | Approximations and complex multiplication according to Ramanujan[END_REF] Strassen's algorithm to the computation of the N -th term of any holonomic sequence in arithmetic complexity quasi-linear in √ N . Our main contribution in this article consists in transferring these results to the q-holonomic framework. It turns out that the resulting algorithms are actually simpler in the q-holonomic case than in the usual holonomic setting, essentially because multipoint evaluation on arithmetic progressions used as a subroutine in Strassen's and Chudnovskys' algorithms is replaced by multipoint evaluation on geometric progressions, which is considerably simpler [START_REF] Bostan | Polynomial evaluation and interpolation on special sets of points[END_REF].

A consequence of our results is that the following apparently unrelated polynomials and rational functions can be evaluated fast (note the change in notation, with the variable x denoted now by q):

• A n (q), the generating function of the number of partitions into n positive integers each occurring at most twice [START_REF] Yang | On the product n≥1 (1 + q n x + q 2n x 2 )[END_REF], i.e., the coefficient of t n in the product k≥1

(1 + q k t + q 2k t 2 ).

• B n (q) := ∞ i=1 (1 -q i ) mod q n ; by Euler's pentagonal theorem [113, §5],

B n (q) = 1 + i(3i+1)<2n (-1) i q i(3 i-1) 2 + q i(3 i+1) 2
.

• The number C n (q) of 2n × 2n upper-triangular matrices over F q (the finite field with q elements), whose square is the zero matrix [START_REF] Kirillov | On a remarkable sequence of polynomials[END_REF]; by [START_REF] Ekhad | The number of solutions of X 2 = 0 in triangular matrices over GF(q)[END_REF], C n (q) is equal to

C n (q) = j 2n n -3j - 2n n -3j -1 • q n 2 -3j 2 -j .
The common feature, exploited by the new algorithm, is that the sequences (A n (q)) n≥0 , (B n (q)) n≥0 , (C n (q)) n≥0 are all q-holonomic. Actually, q-holonomic sequences are ubiquitous, so the range of application of our results is quite broad. This stems from the fact that they are coefficient sequences of power series satisfying q-difference equations, or equivalently, q-shift (or, q-differential) equations. From that perspective, our topic becomes intimately connected with q-calculus. The roots of q-calculus are in works of famous mathematicians such as Rothe [START_REF] Rothe | Formulae de serierum reversione demonstratio universalis signis localibus combinatorico-analyticorum vicariis exhibita[END_REF], Gauss [START_REF] Gauss | Summatio quarundam serierum singularium[END_REF] and Heine [START_REF] Heine | Untersuchungen über die Reihe 1+ (1-q α )(1-q β ) (1-q)(1-q γ ) x+ (1-q α )(1-q α+1 )(1-q β )-q β+1 ) -q)-q 2 )-q γ )(-q γ+1 ) x[END_REF]. The topic gained renewed interest in the first half of the 20th century, with the work, both on the formal and analytic aspects, of Tanner [START_REF] Tanner | On the Enumeration of Groups of Totitives[END_REF], Jackson [START_REF] Jackson | On q-functions and a certain difference operator[END_REF][START_REF] Jackson | On q-definite integrals[END_REF][START_REF] Jackson | q-Difference Equations[END_REF], Carmichael [START_REF] Carmichael | The General Theory of Linear q-Difference Equations[END_REF], Mason [START_REF] Mason | On Properties of the Solutions of Linear q-Difference Equations with Entire Function Coefficients[END_REF], Adams [START_REF] Adams | On the linear ordinary q-difference equation[END_REF][START_REF] Adams | Linear q-difference equations[END_REF], Trjitzinsky [START_REF] Trjitzinsky | Analytic theory of linear q-difference equations[END_REF], Le Caine [START_REF] Caine | The linear q-difference equation of the second order[END_REF] and Hahn [START_REF] Hahn | Über die höheren Heineschen Reihen und eine einheitliche Theorie der sogenannten speziellen Funktionen[END_REF], to name just a few. Modern accounts of the various aspects of the theory (including historical ones) can be found in [START_REF] Vizio | Équations aux q-différences[END_REF][START_REF] Ernst | A comprehensive treatment of q-calculus[END_REF][START_REF] Koekoek | Hypergeometric orthogonal polynomials and their q-analogues[END_REF].

One of the reasons for interest in q-differential equations is that, formally, as q tends to 1, the q-derivative f (qx)-f (x) (q-1)x tends to f ′ (x), thus to every differential equation corresponds a q-differential equation which goes formally to the differential equation as q → 1.

In nice cases, (some of) the solutions of the q-difference equation go to solutions of the associated differential equation as q → 1. An early example of such a good deformation behavior is given by the basic hypergeometric equation of Heine [START_REF] Heine | Untersuchungen über die Reihe 1+ (1-q α )(1-q β ) (1-q)(1-q γ ) x+ (1-q α )(1-q α+1 )(1-q β )-q β+1 ) -q)-q 2 )-q γ )(-q γ+1 ) x[END_REF], see also [95, §1.10].

In computer algebra, q-holonomic sequences were considered starting from the early nineties, in the context of computer-generated proofs of identities in the seminal paper by Wilf and Zeilberger [START_REF] Wilf | An algorithmic proof theory for hypergeometric (ordinary & q) multisum/integral identities[END_REF], notably in Section 5 ("Generalization to q-sums and qmultisums") and in Section 6.4 ("q-sums and integrals"). Creative telescoping algorithms for (proper) q-hypergeometric sequences are discussed in various references [START_REF] Böing | Algorithms for q-hypergeometric summation in computer algebra[END_REF][START_REF] Cartier | Démonstration "automatique" d'identités et fonctions hypergéométriques (d'après D. Zeilberger)[END_REF][START_REF] Petkovšek | [END_REF]; several implementations of those algorithms are described for instance in [START_REF] Kauers | A Mathematica package for q-holonomic sequences and power series[END_REF][START_REF] Paule | A Mathematica q-analogue of Zeilberger's algorithm based on an algebraically motivated approach to q-hypergeometric telescoping[END_REF][START_REF] Riese | qMultiSum-a package for proving q-hypergeometric multiple summation identities[END_REF][START_REF] Sprenger | Algorithmic determination of q-power series for q-holonomic functions[END_REF]. Algorithms for computing polynomial, rational and q-hypergeometric solutions of q-difference equations were designed by Abramov and collaborators [START_REF] Abramov | On polynomial solutions of linear operator equations[END_REF][START_REF] Abramov | Rational solutions of linear difference and q-difference equations with polynomial coefficients[END_REF][START_REF] Abramov | q-hypergeometric solutions of q-difference equations[END_REF][START_REF] Khmel | Improved algorithms for solving difference and q-difference equations[END_REF]. These algorithms are important for several reasons. One is that they lie at the heart of the vast generalization by Chyzak [START_REF] Chyzak | Gröbner bases, symbolic summation and symbolic integration[END_REF][START_REF] Chyzak | An extension of Zeilberger's fast algorithm to general holonomic functions[END_REF] of the Wilf and Zeilberger algorithmic theory, for the treatment of general q-holonomic (not only q-hypergeometric) symbolic summation and integration via creative telescoping. In that context, a multivariate notion of q-holonomy is needed; the foundations of the theory were laid by Zeilberger [START_REF] Zeilberger | A holonomic systems approach to special functions identities[END_REF] and Sabbah [START_REF] Sabbah | Systèmes holonomes d'équations aux q-différences[END_REF] (in the language of D-modules), see also [48, § 2.5] and [START_REF] Garoufalidis | A survey of q-holonomic functions[END_REF].

The simplest non-trivial holonomic sequence is (n!) n≥0 , whose n-th term combinatorially counts the number of permutations of n objects. If instead of direct counting, one assigns to every permutation π its number of inversions inv(π), i.e., the number of pairs 1 ≤ i < j ≤ n with π(i) > π(j), the refined count (by size and number of inversions) is

[n] q ! := (1 + q)(1 + q + q 2 ) • • • (1 + q + • • • + q n-1 ).
This is the q-analogue of n!; it is the simplest non-trivial q-holonomic sequence.

There is also a natural q-analog of the binomial coefficients, called the Gaussian coefficients, defined by

n k q := [n] q ! [k] q ![n -k] q ! .
They have many counting interpretations, e.g., they count the k-dimensional subspaces of F n q (points on Grassmannians over F q ). There are q-analogs to (almost) everything. To select just two more basic examples, the q-analog [14, Thm. 3.3] of the binomial theorem is given by n k=1

(1 + q k-1 x) = n k=0 n k q q ( k 2 ) x k , (1)
and the q-version [14, Thm. 3.4] of the Chu-Vandermonde identity is

n k=0 q k 2 m k q n k q = m + n n q . ( 2
)
The ubiquity of q-holonomic sequences is manifest in plenty of fields: partition theory [START_REF] Andrews | The theory of partitions[END_REF][START_REF] Andrews | The fifth and seventh order mock theta functions[END_REF][START_REF] Liu | Some finite generalizations of Euler's pentagonal number theorem[END_REF][START_REF] Pak | Partition bijections, a survey[END_REF][START_REF] Shanks | A short proof of an identity of Euler[END_REF][START_REF] Yang | On the product n≥1 (1 + q n x + q 2n x 2 )[END_REF] and other subfields of combinatorics [START_REF] Andrews | q-Catalan identities[END_REF][START_REF] Bousquet-Mélou | Convex polyominoes and algebraic languages[END_REF][START_REF] Ekhad | The number of solutions of X 2 = 0 in triangular matrices over GF(q)[END_REF][START_REF] Fürlinger | q-Catalan numbers[END_REF][START_REF] Kirillov | On the number of solutions of the equation X 2 = 0 in triangular matrices over a finite field[END_REF][START_REF] Kirillov | On a remarkable sequence of polynomials[END_REF][START_REF] Yip | Rook placements and Jordan forms of upper-triangular nilpotent matrices[END_REF]; theta functions and modular forms [START_REF] Bellman | A brief introduction to theta functions[END_REF][START_REF] Garvan | New fifth and seventh order mock theta function identities[END_REF][START_REF] Labrande | Computing Jacobi's theta in quasi-linear time[END_REF][START_REF] Labrande | Computing theta functions in quasi-linear time in genus two and above[END_REF][START_REF] Zagier | Elliptic modular forms and their applications[END_REF]; special functions [START_REF] Borwein | Padé approximants for the q-elementary functions[END_REF][START_REF] Ismail | Lectures on q-orthogonal polynomials[END_REF][START_REF] Koekoek | Hypergeometric orthogonal polynomials and their q-analogues[END_REF] and in particular orthogonal polynomials [START_REF] Koornwinder | On Zeilberger's algorithm and its q-analogue[END_REF]; algebraic geometry [START_REF] Ekedahl | Cycle classes on the moduli of K3 surfaces in positive characteristic[END_REF], representation theory [START_REF] Hua | Counting representations of quivers over finite fields[END_REF]; knot theory [START_REF] Garoufalidis | Twisting q-holonomic sequences by complex roots of unity[END_REF][START_REF] Garoufalidis | Irreducibility of q-difference operators and the knot 7 4[END_REF][START_REF] Garoufalidis | The colored HOMFLYPT function is q-holonomic[END_REF][START_REF] Garoufalidis | The colored Jones function is q-holonomic[END_REF][START_REF] Garoufalidis | A survey of q-holonomic functions[END_REF]; Galois theory [START_REF] Hendriks | An algorithm for computing a standard form for second-order linear q-difference equations[END_REF]; number theory [START_REF] Adamczewski | Congruences modulo cyclotomic polynomials and algebraic independence for q-series[END_REF][START_REF] Vizio | Arithmetic theory of q-difference equations: the q-analogue of Grothendieck-Katz's conjecture on p-curvatures[END_REF][START_REF] Osgood | On the diophantine approximation of values of functions satisfying certain linear q-difference equations[END_REF].

The main messages of this article are that for any example of a q-holonomic sequence occurring in those various fields, one can compute selected coefficients faster than by a direct algorithm and that this fact finds a tremendous number of applications.

Complexity basics. We estimate the arithmetic complexities of algorithms by counting arithmetic operations (+, -, ×, ÷) in the base field K at unit cost. We use standard complexity notation, such as M(d) for the cost of degree-d multiplication in K[x], and θ for feasible exponents of matrix multiplication. The best currently known upper bound is θ < 2.3729 [START_REF] Alman | A refined laser method and faster matrix multiplication[END_REF][START_REF] Gall | Powers of tensors and fast matrix multiplication[END_REF]. As usual, O(•) stands for the big-Oh notation and Õ(•) is used to hide polylogarithmic factors in the argument. Most arithmetic operations on univariate polynomials of degree d in K[x] can be performed in quasi-linear complexity Õ(d): multiplication, shift, interpolation, gcd, resultant, etc. A key feature of these results is the reduction to fast polynomial multiplication, which can be performed in time M(d) = O(d log d log log d) [START_REF] Cantor | On fast multiplication of polynomials over arbitrary algebras[END_REF][START_REF] Schönhage | Schnelle Multiplikation von Polynomen über Körpern der Charakteristik 2[END_REF]. Finally, the arithmetic cost of multiplication of polynomial matrices of size n and degree d is denoted by MM(n, d) and we have [START_REF] Bostan | Polynomial evaluation and interpolation on special sets of points[END_REF]. An excellent general reference for these questions is the book by von zur Gathen and Gerhard [START_REF] Zur Gathen | Modern computer algebra[END_REF].

MM(n, d) = O(n θ d + n 2 M(d)) = Õ(n θ d)
A short version of this article has appeared at the ISSAC'20 conference [START_REF] Bostan | Computing the N -th Term of a q-Holonomic Sequence[END_REF]. In the present version, we included the proofs of Theorems 6 and 8, we added a new Theorem 5 containing a detailed complexity analysis of the main algorithm (Algorithm 3) with respect to all parameters, and we displayed pseudo-code for the algorithms as well as figures visualizing their performance. We also elaborated on a task which was mentioned as future work in the previous version, namely the application of our methods to the computation of curvatures of q-difference equations, see §4.4.

The structure of the article is as follows: in Section 2 we deal with the tasks of evaluating Q N (x) and R N (x). We show that these are two instances of the same problem and provide Algorithm 3 which solves both in O(M( √ N )) arithmetic complexity. Section 3 is devoted to the main results; we prove there that Algorithm 3 can be used for computing terms of any q-holonomic sequence with the same cost, and provide extensions and more insight. In the same section we also consider the bit-complexity model. We identify and elaborate on several applications for our result in Section 4. In Section 5 we report on implementations of our algorithms, which deliver encouraging timings, and we finally describe future tasks and investigation fields in Section 6.

Two motivating examples

Before presenting our main results in Section 3, we describe in this section the approach and main ideas on two basic examples. Both examples concern the fast evaluation of special families of univariate polynomials. In §2.1, we consider polynomials of the form ℓ (x -q ℓ ), and in §2.2 sparse polynomials of the form ℓ p ℓ x aℓ 2 +bℓ . In both cases, we first present fast ad-hoc algorithms, then introduce equally fast alternative algorithms, which have the nice feature that they will be generalizable to a broader setting.

Evaluation of some structured polynomials

Here is our first example, that emerged from a question asked to the first author by Luca De Feo (private email communication, 10 January 2020); this was the starting point of the article.

Let q be an element of the field K, and consider the polynomial

F (x) := N -1 i=0 (x -q i ) ∈ K[x]. (3)
Given another element α ∈ K, how fast can one evaluate F (α)?

If q = 0, then F (α) = α N can be computed in O(log N ) operations in K, by binary powering. We assume in what follows that q is nonzero. Obviously, a direct algorithm consists in computing the successive powers q, q 2 , . . . , q N -1 using O(N ) operations in K, then computing the elements α -1, α -q, . . . , α -q N -1 in O(N ) more operations in K, and finally returning their product. The total arithmetic cost of this algorithm 1 is O(N ), linear in the degree of F . Is it possible to do better? The answer is positive, as one can use the following babystep/giant-step strategy, in which, in order to simplify things, we assume that N is a perfect square, N = s 2 .

Algorithm 1

1. (Baby-step) Compute the values of q, q 2 , . . . , q s-1 , and deduce the coefficients of the polynomial

G(x) := s-1 j=0 (x -q j ).
2. (Giant-step) Compute Q := q s , Q 2 , . . . , Q s-1 , and deduce the coefficients of the polynomial

H(x) := s-1 k=0 (α -Q k • x).

Return the resultant Res(G, H).

By the basic property of resultants, the output of this algorithm is

Res(G, H) = s-1 j=0 H(q j ) = s-1 j=0 s-1 k=0 α -q sk+j = N -1 i=0 (α -q i ) = F (α).
Using the fast subproduct tree algorithm [START_REF] Zur Gathen | Modern computer algebra[END_REF]Algorithm 10.3], one can perform the baby-step (1) as well as the giant-step [START_REF] Abramov | Rational solutions of linear difference and q-difference equations with polynomial coefficients[END_REF] in O(M( √ N ) log N ) operations in K, and by [START_REF] Zur Gathen | Modern computer algebra[END_REF]Corollary 11.19] the same cost can be achieved for the resultant computation in step [START_REF] Abramov | A direct algorithm to compute rational solutions of first order linear q-difference systems[END_REF]. Using fast polynomial multiplication, we conclude that F (α) can be computed in arithmetic complexity quasi-linear in √ N . Note that if N is not a perfect square, then one can compute F (α) as

F (α) = F 1 (α)F 2 (α), where F 1 (α) := ⌊ √ N ⌋ 2 -1 i=0 (α -q i ) is computed as in Algorithm 1, while F 2 (α) := N -1 i=⌊ √ N ⌋ 2 (α -q i ) can be computed naively, since N -⌊ √ N ⌋ 2 = O( √ N ).
It is possible to speed up the previous algorithm by a logarithmic factor in N using a slightly different scheme, still based on a baby-step/giant-step strategy, but exploiting the fact that the roots of F are in geometric progression. Again, we assume that N = s 2 is a perfect square. This alternative algorithm goes as follows. Note that it is very close in spirit to Pollard's algorithm described on page 523 of [START_REF] Pollard | Theorems on factorization and primality testing[END_REF].

Algorithm 2

1. (Baby-step) Compute q, q 2 , . . . , q s-1 , and deduce the coefficients of the polynomial P (x) := s-1 j=0 (α -q j • x). 2. (Giant-step) Compute Q := q s , Q 2 , . . . , Q s-1 , and evaluate P simultaneously at

1, Q, . . . , Q s-1 . 3. Return the product P (Q s-1 ) • • • P (Q)P (1).
Obviously, the output of this algorithm is

s-1 k=0 P (Q k ) = s-1 k=0 s-1 j=0 (α -q j • q sk ) = N -1 i=0 (α -q i ) = F (α).
As pointed out in the remarks after the proof of [42, Lemma 1], one can compute P (x) = P s (x) = s-1 j=0 (α -q j • x) in step (1) without computing the subproduct tree, by using a divide-and-conquer scheme which exploits the fact that P 2t (x) = P t (q t x) • P t (x) and P 2t+1 (x) = (α-q 2t x)•P t (q t x)•P t (x). The cost of this algorithm is O(M( √ N )) operations in K.

As for step [START_REF] Abramov | Rational solutions of linear difference and q-difference equations with polynomial coefficients[END_REF], one can use the fast chirp transform algorithms of Rabiner, Schafer and Rader [START_REF] Rabiner | The chirp z-transform algorithm and its application[END_REF] and of Bluestein [START_REF] Bluestein | A linear filtering approach to the computation of the discrete Fourier transform[END_REF]. These algorithms rely on the following observation: writing

Q ij = Q ( i+j 2 ) •Q -( i 2 ) •Q -( j 2
) and P (x) = s j=0 c j x j implies that the needed values

P (Q i ) = s j=0 c j Q ij , 0 ≤ i < s, are P (Q i ) = Q -( i 2 ) • s j=0 c j Q -( j 2 ) • Q ( i+j 2 ) , 0 ≤ i < s,
in which the sum is simply the coefficient of x s+i in the product

  s j=0 c j Q -( j 2 ) x s-j   2s ℓ=0 Q ( ℓ 2 ) x ℓ .
This polynomial product can be computed in 2 M(s) operations (and even in M(s)+O(s) using the transposition principle [START_REF] Bostan | Tellegen's principle into practice[END_REF][START_REF] Hanrot | The middle product algorithm. I[END_REF], since only the median coefficients x s , . . . , x 2s-1 are actually needed). In conclusion, step (2) can also be performed in O(M( √ N )) operations in K, and thus O(M( √ N )) is the total cost of this second algorithm. We have chosen to detail this second algorithm for several reasons: not only because it is faster by a factor log(N ) compared to the first one, but more importantly because it has a simpler structure, which will be generalizable to the general q-holonomic setting. In fact, we do not provide a pseudo-code implementation for this algorithm, since we will do so for the more general case (Algorithm 3).

Evaluation of some sparse polynomials

Let us now consider the sequence of sparse polynomial sums

v (p,a,b) N (q) = N -1 n=0 p n q an 2 +bn ,
where p ∈ K and a, b ∈ Q such that 2a, a + b are both integers. Typical examples are (truncated) modular forms [START_REF] Paule | Rogers-Ramanujan functions, modular functions, and computer algebra[END_REF], which are ubiquitous in complex analysis [START_REF] Bellman | A brief introduction to theta functions[END_REF], number theory [START_REF] Zagier | Elliptic modular forms and their applications[END_REF] and combinatorics [START_REF] Andrews | The theory of partitions[END_REF]. For instance, the Jacobi theta function ϑ 3 depends on two complex variables z ∈ C, and τ ∈ C with ℑ(τ ) > 0, and it is defined by

ϑ 3 (z; τ ) = ∞ n=-∞ e πi(n 2 τ +2nz) = 1 + 2 ∞ n=1 η n q n 2 ,
where q = e πiτ is the nome (|q| < 1) and η = e 2πiz . Here, K = C. Another example is the Dedekind eta function, appearing in Euler's famous pentagonal theorem [113, §5], which has a similar form

q 1 24 • 1 + ∞ n=1 (-1) n q n(3 n-1) 2 + q n(3 n+1) 2
, with q = e 2πiτ .

Moreover, sums of the form v

(1,a,b) N (q) = N -1 n=0 q an 2 +bn , over K = Q or K = F 2 ,
crucially occur in a recent algorithm by Tao, Crott and Helfgott [START_REF] Tao | Deterministic methods to find primes[END_REF] for the efficient construction of prime numbers in given intervals, e.g., in the context of effective versions of Bertrand's postulate. Actually, (the proof of) Lemma 3.1 in [START_REF] Tao | Deterministic methods to find primes[END_REF] contains the first sublinear complexity result for the evaluation of the sum v (p,a,b) N (q) at an arbitrary point q; namely, the cost is O(N θ/3 ), where θ ∈ [2, 3] is any feasible exponent for matrix multiplication. Subsequently, Nogneng and Schost [START_REF] Nogneng | On the evaluation of some sparse polynomials[END_REF] designed a faster algorithm, and lowered the cost down to Õ( √ N ). Our algorithm is similar in spirit to theirs, as it also relies on a baby-step/giant-step strategy.

Let us first recall the principle of the Nogneng-Schost algorithm [START_REF] Nogneng | On the evaluation of some sparse polynomials[END_REF]. Assume as before that N is a perfect square, N = s 2 . The starting point is the remark that

v (p,a,b) N (q) = N -1 n=0 p n q an 2 +bn = s-1 k=0 s-1 j=0
p j+sk q a(j+sk) 2 +b(j+sk) can be written s-1 k=0 p sk q as 2 k 2 +bsk • P (q 2ask ), where P (y) := s-1 j=0 p j q aj 2 +bj y j .

Therefore, the computation of v (p,a,b) N (q) can be reduced essentially to the simultaneous evaluation of the polynomial P at s = 1 + deg(P ) points (in geometric progression), with arithmetic cost O(M( √ N )). We now describe an alternative algorithm, of similar complexity O(M( √ N )), with a slightly larger constant in the big-Oh estimate, but whose advantage is its potential of generality.

Let us denote by u n (q) the summand p n q an 2 +bn . Clearly, the sequence (u n (q)) n≥0 satisfies the recurrence relation

u n+1 (q) = A(q, q n ) • u n (q), where A(x, y) := px a+b y 2a .
As an immediate consequence, the sequence with general term v n (q) := n-1 k=0 u k (q) satisfies a similar recurrence relation

v n+2 (q) -v n+1 (q) = A(q, q n ) • (v n+1 (q) -v n (q)), (4)
with initial conditions v 0 (q) = 0 and v 1 (q) = 1. This scalar recurrence of order two is equivalent to the first-order matrix recurrence

v n+2 v n+1 = A(q, q n ) + 1 -A(q, q n ) 1 0 × v n+1 v n .
By unrolling this matrix recurrence, we deduce that

v n+1 v n = M (q n-1 ) v n v n-1 = M (q n-1 ) • • • M (q)M (1) × 1 0 ,
where

M (x) := pq a+b x 2a + 1 -pq a+b x 2a 1 0 , hence v N = 0 1 × M (q N -1 ) • • • M (q)M (1) × 1 0
. Therefore, the computation of v N reduces to the computation of the "matrix q-factorial" M (q N -1 ) • • • M (q)M (1), which can be performed fast by using a baby-step/giant-step strategy similar to the one of the second algorithm in §2.1. Again, we assume for simplicity that N = s 2 is a perfect square. The algorithm goes as follows.

Algorithm 3 (matrix q-factorial)

(1) (Baby-step) Compute q, q 2 , . . . , q s-1 ; deduce the coefficients of the polynomial matrix

P (x) := M (q s-1 x) • • • M (qx)M (x). (2) (Giant-step) Compute Q := q s , Q 2 , . . . , Q s-1
, and evaluate (the entries of) P (x)

simultaneously at 1, Q, . . . , Q s-1 . (3) Return the product P (Q s-1 ) • • • P (Q)P (1).
Clearly, this algorithm generalizes Algorithm 2 in §2.1 and, as promised, we also provide a detailed pseudo-code implementation: Step1 and Step2 & Step3:

Algorithm 3 (Step1) Input: s, q, M (x) Output: M (q s-1 x) • • • M (qx)M (x) 1: q s ← [q, q 2 , . . . , q s-1 ] 2: t ← s 3: function BS(t) 4: if t = 1 then 5: return M (x) 6:
if t is even then 7: p 2 (x) ← p 1 (q (t-1)/2 x) ⊲ Using q s 13:

p 1 (x) ← BS(t/2) 8: p 2 (x) ← p 1 (q t/2 x) ⊲ Using q s 9: return p 2 (x) • p 1 (x) ⊲ Fast
p 3 (x) ← M (q t-1 x) ⊲ Using q s 14: return p 3 (x) • p 2 (x) • p 1 (x) ⊲ Fast polynomial multiplication Algorithm 3 (Step2 & Step3) Input: s, Q, P (x) Output: P (Q s-1 ) • • • P (1)
Assumptions: Q = 0, P (x) polynomial matrix of size n × n and degree d ≥ s.

1: Q d ← [Q, Q 2 , . . . , Q d-1 ] 2: Q ′ ← 1/Q 3: Q ′ d ← [Q -( d 2 ) , . . . , Q -( 1 2 ) , Q ( 0 2 ) , . . . , Q ( 2d 2 ) ]
⊲ Using Q d and Q ′ 4: P s-1 , . . . , P 0 ⊲ Empty n × n matrices 5: for i from 1 to n do 6:

for j from 1 to n do 7:

p(x) ← P (x) i,j ⊲ p(x) = c 0 + c 1 x + • • • + c d x d 8: p 1 (x) ← d ℓ=0 c ℓ Q -( ℓ 2 ) x d-ℓ ⊲ Using Q ′ d 9: p 2 (x) ← 2d ℓ=0 Q ( ℓ 2 ) x ℓ ⊲ Using Q ′ d 10: p 3 (x) ← P 1 (x) • P 2 (x) ⊲ p 3 (x) = 3d ℓ=0 r ℓ x ℓ ; fast multiplication 11:
for k from 0 to s -1 do 12: P ← P k • P return P By the same observations as in Algorithm 2 in §2.1, the complexity of Algorithm 3 already is quasi-linear in √ N . In the next section we will discuss the complexity not only with respect to N , but to the matrix size and degree as well.

(P k ) i,j ← r d+k ⊲ P ℓ = P (Q ℓ ) for ℓ = 0, . . . ,
We remark that when applied to the computation of v

M (x) := prx + 1 -prx 1 0 , with r := q a+b , (5)
and if q := q 2a , then the following matrix q-factorials coincide:

M (q N -1 ) • • • M (q)M (1) = M (q N -1 ) • • • M (q) M (1)
.

Algorithm 4 (matrix q-factorial, variant) (0) (Precomputation) Compute r := q a+b , q := q 2a , and M in ( 5).

(1) (Baby-step) Compute q, q2 , . . . , qs-1 ; deduce the coefficients of the polynomial matrix

P (x) := M (q s-1 x) • • • M (qx) M (x).
(2) (Giant-step) Compute Q := qs , Q2 , . . . , Qs-1 , and evaluate (the entries of) P (x) simultaneously at 1, Q, . . . , Qs-1 . (3) Return the product P ( Qs-1 )

• • • P ( Q) P (1).
Using binary powering, the cost of the additional precomputation in step (0) is only logarithmic in a and b. In exchange, the new steps (2) and (3) are performed on matrices whose degrees do not depend on a and b anymore (in the previous, unoptimized, version the degrees of the polynomial matrices were linear in a and b). The total arithmetic cost with respect to N is still quasi-linear in √ N .

In the next section, we will show that Algorithm 3 can be employed for the fast computation of the N -th term of any q-holonomic sequence. Note that the trick in Algorithm 4 relies on the fact that M (x), coming from the recurrence for v (p,a,b) N (q), contains only pure powers of x and q. We cannot hope for this phenomenon in general, however we advise to bear this simplification in mind for some practical purposes. In any case, we can improve on the quasi-linear cost in the degree d of the polynomial matrix M (x) in Algorithm 3, obtaining a complexity of essentially √ d; in essence, the idea consists in choosing s = N/d rather than √ N , see §3.4.

Main results

In this section, we generalize the algorithms from §2, and show that they apply to the general setting of q-holonomic sequences.

Preliminaries

A sequence u n = u n (q) is q-holonomic if it satisfies a nontrivial q-recurrence, that is, a linear recurrence with coefficients given by polynomials in q and q n . Definition 1 (q-holonomic sequence). Let K be a field, and q ∈ K. A sequence (u n (q)) n≥0 in K N is called q-holonomic if there exist r ∈ N and polynomials c 0 (x, y), . . . , c r (x, y) in K[x, y], with c r (x, y) = 0, such that c r (q, q n )u n+r (q) + • • • + c 0 (q, q n )u n (q) = 0, for all n ≥ 0.

(

) 6
The integer r is called the order of the q-recurrence (6). When r = 1, we say that (u n (q)) n≥0 is q-hypergeometric.

The most basic examples are the q-bracket and the q-factorial,

[n] q := 1 + q + • • • + q n-1 and [n] q ! := n k=1 [k] q . ( 7
)
They are clearly q-holonomic, and even q-hypergeometric.

The sequences

(u n ) n≥0 = (q n ) n≥0 , (v n ) n≥0 = (q n 2 ) n≥0 and (w n ) n≥0 = (q ( n 2
) ) n≥0 are also q-hypergeometric, since they satisfy the recurrence relations

u n+1 -qu n = 0, v n+1 -q 2n+1 v n = 0, w n+1 -q n w n = 0.
However, the sequence (q n 3 ) n≥0 is not q-holonomic [72, Ex. 2.2(b)]. More generally, this also holds for the sequence (q n s ) n≥0 , for any s > 2, see [START_REF] Bézivin | Sur les équations fonctionnelles aux q-différences[END_REF]Th. 4.1] and also [START_REF] Garoufalidis | The degree of a q-holonomic sequence is a quadratic quasi-polynomial[END_REF]Th. 1.1].

Another basic example is the q-Pochhammer symbol

(x; q) n := n-1 k=0 (1 -xq k ), (8)
which is also q-hypergeometric, since (x; q) n+1 -(1 -xq n )(x; q) n = 0. In particular, the sequence (q; q) n := n k=1 (1 -q k ), also denoted (q) n , is q-hypergeometric and satisfies (q) n+1 -(1 -q n+1 )(q) n = 0. In Section §2 we encountered v (p,a,b) n = n k=0 p k q ak 2 +bk , which is q-holonomic (see Eq. ( 4)), but generally not q-hypergeometric.

Note that (6) reduces to a C-linear recurrence, i.e. a linear recurrence with constant coefficients, if all polynomials c 0 (x, y), . . . , c r (x, y) are constant in the variable y. For these kinds of sequences there exist quasi-optimal algorithms [START_REF] Bostan | A Simple and Fast Algorithm for Computing the N -th Term of a Linearly Recurrent Sequence[END_REF][START_REF] Fiduccia | An efficient formula for linear recurrences[END_REF][START_REF] Miller | An algorithm for evaluation of remote terms in a linear recurrence sequence[END_REF], therefore we assume from now on that the maximal degree d of c 0 (x, y), . . . , c r (x, y) in y is positive.

As mentioned in the introduction, q-holonomic sequences show up in various contexts. As an example, in (quantum) knot theory, the ("colored") Jones function of a (framed oriented) knot (in 3-space) is a powerful knot invariant, related to the Alexander polynomial [START_REF] Bar-Natan | On the Melvin-Morton-Rozansky conjecture[END_REF]; it is a q-holonomic sequence of Laurent polynomials [START_REF] Garoufalidis | The colored Jones function is q-holonomic[END_REF]. Its recurrence equations are themselves of interest, as they are closely related to the A-polynomial of a knot, via the AJ conjecture [START_REF] Detcherry | A diagrammatic approach to the AJ conjecture[END_REF][START_REF] Garoufalidis | On the characteristic and deformation varieties of a knot[END_REF][START_REF] Garoufalidis | Quantum knot invariants[END_REF], verified in some cases using massive computer algebra calculations [START_REF] Garoufalidis | Irreducibility of q-difference operators and the knot 7 4[END_REF].

It is well known that the class of q-holonomic sequences is closed under several operations, such as addition, multiplication, Hadamard product and monomial substitution [START_REF] Garoufalidis | A survey of q-holonomic functions[END_REF][START_REF] Kauers | A Mathematica package for q-holonomic sequences and power series[END_REF][START_REF] Koepf | Properties of q-holonomic functions[END_REF]. All these closure properties are effective, i.e., they can be executed algorithmically on the level of q-recurrences. Several computer algebra packages are available for the manipulation of q-holonomic sequences, e.g., the Mathematica packages qGeneratingFunctions [START_REF] Kauers | A Mathematica package for q-holonomic sequences and power series[END_REF] and HolonomicFunctions [START_REF] Koutschan | A fast approach to creative telescoping[END_REF], and the Maple packages qsum [START_REF] Böing | Algorithms for q-hypergeometric summation in computer algebra[END_REF], qFPS [START_REF] Sprenger | Algorithmic determination of q-power series for q-holonomic functions[END_REF], qseries and QDifferenceEquations.

A simple but useful fact is that the order-r scalar q-recurrence (6) can be translated into a first-order recurrence on r × 1 vectors:

   u n+r . . . u n+1    =      -cr-1 cr • • • -c1 cr -c0 cr 1 • • • 0 0 . . . . . . . . . . . . 0 • • • 1 0      ×    u n+r-1 . . . u n    . (9)
In particular, the N -th term of the q-holonomic sequence (u n ) is simply expressible in terms of the matrix q-factorial

M (q N -1 ) • • • M (q)M (1), (10)
where M (q n ) denotes the companion matrix from equation ( 9). This observation is crucial, since it exposes the connection to the algorithms presented in the previous section.

Computation of the q-factorial

We now give the promised q-analogue of Strassen's result on the computation of N ! in O(M( √ N ) log N ) arithmetic operations. Note that Strassen's case q = 1 is also covered by [38, §6], where the cost O(M( √ N )) is reached under some invertibility assumptions.

Theorem 1. Let K be a field, let q ∈ K \ {1} and N ∈ N. The q-factorial [N ] q ! can be computed using O(M( √ N )) operations in K. The same is true for the q-Pochhammer symbol (α; q) N for any α ∈ K.

Proof. If α = 0, then (α; q) N = 1. If q = 0, then [N ] q ! = 1 and (α; q) N = 1 -α. We can assume that q ∈ K \ {0, 1} and α ∈ K \ {0}. We have [N ] q ! = r N • F (q -1 ) and (α; q) N = α N • F (α -1 ), where r := q/(1 -q) and F (x) := N -1 i=0 (x -q i ). Algorithm 2 can be used to compute F (q -1 ) and F (α -1 ) in O(M( √ N )) operations in K. The cost of computing r N and α N is O(log N ), and thus it is negligible.

Corollary 2. Under the assumptions of Theorem 1 and for any n ∈ N, one can compute in O(M( √ N )) operations in K:

• the q-binomial coefficient N n q ;

• the coefficient of x n in the polynomial

N k=1 (1 + q k-1 x);
• the sum N -n 0 q n 0 q + q N -n

1 q n 1 q + • • • + q n 2 N -n n q n n q .
Proof. The first assertion is a direct consequence of Theorem 1. The second assertion is a consequence of the first one, and of (1). The third assertion is a consequence of the first one, and of (2).

3.3. N -th term of a q-holonomic sequence

We now offer the promised q-analogue of Chudnovskys' result on the computation of the N -th term of an arbitrary holonomic sequence in O(M( √ N ) log N ) arithmetic operations. Note that Chudnovskys' case q = 1 is also covered by [38, §6], where the improved cost O(M( √ N )) is reached under additional invertibility assumptions.

Theorem 3. Let K be a field, q ∈ K \ {1} and N ∈ N. Let (u n (q)) n≥0 be a q-holonomic sequence satisfying recurrence [START_REF] Adamczewski | Congruences modulo cyclotomic polynomials and algebraic independence for q-series[END_REF], and assume that c r (q, q k ) is nonzero for k = 0, . . . , N -1. Then, u N (q) can be computed in O(M( √ N )) operations in K.

Proof. Using equation ( 9), it is enough to show that the matrix q-factorial

M (q N -1 ) • • • M (q)M (1) can be computed in O(M( √ N ))
, where M (q n ) denotes the companion matrix from equation [START_REF] Al-Salam | Some orthogonal q-polynomials[END_REF]. Algorithm 3 adapts mutatis mutandis to this effect.

Remark that if q is a root of unity of order n < N , then the computation of U N (q) = M (q N -1 ) • • • M (q)M (1) can be simplified using

U N (q) = M (q k ) • • • M (1) • U n (q) r ,
where r = ⌊(N -1)/n⌋ and k = N -1 -rn. Algorithm 3 is used to compute U n (q) and then its r-th power is then deduced via binary powering. Finally, the product M (q k ) • • • M (1) is again computed using Algorithm 3. The total cost therefore consists of just O(M( √ n) + log(N )) arithmetic operations. It follows that if, for instance, the base field K is the prime field F p , then the prime number p should be larger than N in order to exhibit the full strength of the presented algorithms.

Corollary 4. Let K be a field, q ∈ K not a root of unity, and N ∈ N. Let e q (x) be the q-exponential series e q (x) := n≥0

x n

[n] q ! , and let E (N ) q (x) := e q (x) mod x N be its polynomial truncation of degree N -1.

If α ∈ K, then one can compute E (N ) q (α) in O(M( √ N )) operations in K.
Proof. Denote the summand α n

[n]q! by u n (q). Then (u n (q)) n≥0 is q-hypergeometric, and satisfies the recurrence [n+1] q u n+1 (q)-αu n = 0, therefore v N (q) := N -1 i=0 u i (q) satisfies the second-order recurrence [n+1] q (v n+2 (q)-v n+1 (q))-α(v n+1 (q)-v n (q)) = 0. Applying Theorem 3 to v N (q) concludes the proof.

The same result holds true if e q (x) is replaced by any power series satisfying a qdifference equation. For instance, one can evaluate fast all truncations of Heine's qhypergeometric series

2 φ 1 ([a, b], [c]; q; x) := n≥0 (a; q) n (b; q) n (c; q) n •
x n (q) n .

Complexity analysis and computation of several terms

Theorem 3 established an O(M( √ N )) cost of the presented method for computing the N -th term of a q-holonomic sequence. We now aim at performing a detailed complexity analysis with respect to all input parameters. So we need to discuss the complexity of Algorithm 3, where we assume that M (x) ∈ M n (K[x] d ) is an n × n polynomial matrix of degree d ≥ 1. We wish to examine the amount of field operations in K needed for the computation of M (q N -1 ) • • • M (1) in terms of N, d and n. Recall that MM(n, d) controls the arithmetic complexity of the product in

M n (K[x] d ) and it holds that MM(n, d) = O(n θ d + n 2 M(d)) = Õ(n θ d).
First, we will examine the direct application of Algorithm 3, where s = √ N , to M (x). It turns out that the dominating part is step [START_REF] Abramov | On polynomial solutions of linear operator equations[END_REF], where we compute P s (x) = M (q s-1 x) • • • M (x) using the divide-and-conquer scheme P 2t (x) = P t (q t x) • P t (x) and P 2t+1 (x) = M (q 2t x) • P t (q t x) • P t (x). Note that P t (x) is an n × n polynomial matrix of degree at most td and therefore the cost of this step is

O(MM(n, sd)) = Õ(n θ d √ N ).
Step ( 2) is done component-wisely at each entry of P (x). By the explained fast chirp transform algorithms, it essentially boils down to n 2 multiplications of two polynomials, one of degree sd and the other of degree 2sd. The cost of the second step is therefore

O(n 2 M(sd)) = Õ(n 2 d √ N )
. The last step is the multiplication of N/s = s matrices with entries in K and has therefore an arithmetic complexity of O(n θ s) = Õ(n θ √ N ). If d < N is a parameter of interest, then there is a better choice of s rather than √ N . We saw that the polynomial P s (x) has degree sd and we must evaluate it at N/s points. The optimal pick for s is therefore s = N/d, which we again can assume to be integer 2 . Then, by the same arguments as above, the costs of the three steps are O(MM(n,

√ N d)) = Õ(n θ √ N d), O(n 2 M( √ N d)) = Õ(n 2 √ N d) and O(n θ √ N d) respec- tively.
Now, we address specifically the computation of the N -th term in a q-holonomic sequence. If (u n (q)) n≥0 is given by a q-recurrence c r (q, q n )u n+r (q) + • • • + c 0 (q, q n )u n (q) = 0, for q ∈ K and for polynomials c j (x, y) ∈ K[x, y], then as observed before, we can compute u N (q) via 1 c r (q, q N -1 ) • • • c r (q, q)c r (q, 1)

• 0 • • • 0 1 × M (q N -1 ) • • • M (q) M (1) ×    1 . . . 0    ,
where now

M (x) := c r (q, x) • M (x) =      -c r-1 (q, x) • • • -c 1 (q, x) -c 0 (q, x) c r (q, x) • • • 0 0 . . . . . . 0 • • • c r (q, x) 0     
.

2 Similarly as before, if N/d is not an integer, then we can compute u N 1 (q) first, where N 1 = ⌊ N/d⌋ 2 d, and then proceed "naively". Note that N -

N 1 < 2 √ N d -d = O( √ N d).
Hence, we are interested in M (q N -1 ) • • • M (1) and c r (q, q N -1 ) • • • c r (q, 1). If the degrees of c 0 (q, y), . . . , c r (q, y) are bounded by d, then the considerations above imply that the two q-factorials can be computed in O(MM(r,

√ N d)+r 2 M( √ N d))
and O(M( √ N d)) operations in K, respectively. We obtain the following theorem (compare with [37, Thm. 2]).

Theorem 5. Under the assumptions of Theorem 3, let d ≥ 1 be the maximum of the degrees of c 0 (q, y), . . . , c r (q, y). Then, for any N > d, the term u N (q) can be computed in

O(r θ √ N d + r 2 M( √ N d)) operations in K.
Theorem 3 can be adapted to the computation of several coefficients of a q-holonomic sequence. The proof is similar to that of Theorem 15 in [START_REF] Bostan | Linear recurrences with polynomial coefficients and application to integer factorization and Cartier-Manin operator[END_REF], however simpler, because we deal with geometric progressions instead of arithmetic ones. Theorem 6. Under the assumptions of Theorem 5, let

N 1 < N 2 < • • • < N n = N be positive integers, where n ≤ √ N .
Then, the terms u N1 (q), . . . , u Nn (q) can be computed altogether in O(M( √ N ) log N ) operations in K.

Proof. As before, we assume that N is a perfect square; let s = √ N . Examining the presented algorithms, we notice that on the way of computing the matrix q-factorial U N := M (q N -1 ) • • • M (q)M (1) we obtain the evaluated polynomials P (1), P (Q), . . . , P (Q s-1 ), where Q := q s and P (x) := M (q s-1 x) • • • M (qx)M (x). The q-factorial U N is then found by step (3) by trivially multiplying P (Q s-1 ) • • • P (Q)P [START_REF] Abramov | On polynomial solutions of linear operator equations[END_REF]. Observe that while multiplying together from right to left we actually also automatically compute

P (Q j-1 ) • • • P (Q)P (1) = M (q sj-1 ) • • • M (q)M (1) = U sj ,
for every j = 1, . . . , s. It follows that employing Algorithm 3 and by simply taking the top right element of each U sj , we find not only u N , but actually u s , u 2s , . . . , u s 2 = u N . This already indicates that simultaneous computation of s terms is achievable in similar complexity after some "distillation". In general, we are interested in the sequence of u i (q) at indices i = N 1 , . . . , N n , hence we need to perform the following refinement step.

Let d 0 ∈ N be a positive integer with d 0 ≤ n ≤ √ N and assume that for some k

(0) 1 , . . . , k (0) n with k (0) j ≤ N j < k (0) j + 2d 0 we already know the values U k (0) 1 , . . . , U k (0) n .
Then we can use a similar strategy as in step (1) of Algorithm 3 and deduce the polynomial matrix

P d0 (x) = M (q d0-1 x) • • • M (qx)M (x). Compute then the values q k (0)
1 , . . . , q k (0) n and evaluate P d0 (x) simultaneously at them. For each j = 1, . . . , n it holds that

P d0 (q k (0) j ) • U k (0) j = U k (0) j +d0 .
We perform this multiplication for those indices j for which k

(0) j + d 0 ≤ N j < k (0)
j + 2d 0 . For these j we then set k

O(M( √ N )) + O(M( √ N ) log N ) = O(M( √ N ) log N ).
The same idea applies in the following corollary which states that if n < √ N /N ε for some ε > 0, then we can omit the log-factor in N : Corollary 7. Under the assumptions of Theorem 5, let

N 1 < N 2 < • • • < N n = N be pos- itive integers, where n < N 1 2 -ε for some 0 < ε < 1 2 .
Then, the terms u N1 (q), . . . , u Nn (q) can be computed altogether in O(M( √ N )) operations in K.

Proof. Here, we follow the exact same procedure as in the proof before. Then, regarding complexity, we use

O(M(n)) = O(M(N 1 2 -ε )) = O(M( √ N )N -ε
) and obtain that the same method yields a total arithmetic cost of O(M(

√ N )) + O(M( √ N )N -ε log N ) = O(M( √ N )).
Remark that regarding a detailed complexity analysis for the computation of several coefficients, we have the following trade-off: either we compute at most √ N terms in the arithmetic complexity O((r

θ d √ N + r 2 M(d √ N )) log N ), or at most N/d terms, but in a better cost of O((r θ √ N d + r 2 M( √ N d)) log N ).
The proofs combine the considerations above with setting s = √ N and s = N/d respectively. In both cases we can get rid of the log-factor in N like in Corollary 7 by computing a factor of N ε less terms.

The case q is an integer: bit complexity

Until now, we only considered the arithmetic complexity model, which is very wellsuited to measure the algorithmic cost when working in algebraic structures whose basic internal operations have constant cost (such as finite fields, or floating point numbers). Now we discuss here the case where q is an integer (or rational) number. The arithmetic complexity model needs to be replaced by the bit-complexity model.

Recall that the most basic operations on integer numbers can be performed in quasioptimal time, that is, in a number of bit operations which is almost linear, up to logarithmic factors, in (the maximum of) their bit size. The most basic operation is integer multiplication, for which quasi-linear time algorithms are known since the early seventies, starting with the famous paper by Schönhage and Strassen [START_REF] Schönhage | Schnelle Multiplikation grosser Zahlen[END_REF] who showed that two n-bit integers can be multiplied in O(n log n log log n) bit operations. After several successive improvements, e.g., [START_REF] Fürer | Faster integer multiplication[END_REF][START_REF] Harvey | Even faster integer multiplication[END_REF], we know as of 2020 that two n-bit integers can be multiplied in time O(n log n) [START_REF] Harvey | Integer multiplication in time O(n log n)[END_REF]. We shall call the cost of multiplying two n-bit integers M Z (n).

In this context, the matrix q-factorials from §3.1 are computed by binary splitting rather than by baby-step/giant-steps. Recall that this phenomenon already occurs in the usual holonomic setting. For example, the bitsize of u N = N ! is O(N log N ), however both, the "naive" method of computing it using u n = nu n-1 , or Strassen's babysteps/giant-steps method, yield worse bit complexity. In the naive approach the problem is that integers of unbalanced bitsize are multiplied together and hence not the full power of fast integer multiplication techniques can be employed. A more clever and very simple way is to just use the fact that

N ! = (1 • • • ⌊N/2⌋) × ((⌊N/2⌋ + 1) • • • N ),
and that the bitsizes of both factors have magnitude O(N/2 log(N )) = Õ(N ). Thus, fast integer multiplication can be used to multiply them in Õ(N ) bit-complexity. This idea results in the binary splitting Algorithm 5. This algorithm is very classical, and we only recall it for completeness. See [23, §12] for a good survey on this technique, and its applications.

Algorithm 5 (BinSplit)

Input: A = [a 1 , . . . , a N ]
list of elements from some arbitrary ring R Output:

a N • • • a 1 1: function F (A) 2: if N = 1 then 3: return A[1] 4: return F (A[⌊N/2⌋ + 1, . . . , N ]) • F (A[1, . . . , ⌊N/2⌋])
Assume that each a i in the input of BinSplit has at most k bits and let C(n) be the complexity of

BinSplit if A = [a 1 , . . . , a n ] is n-dimensional. It follows that C(N ) ≤ 2C(⌈N/2⌉) + M Z (N/2 • k),
where M Z (n) = Õ(n) is the cost of multiplication of integers with n bits. We obtain C(N ) = Õ(N k). Hence, using this method, the computation of u N = N ! has Õ(N ) bit-complexity, which is quasi-optimal. Moreover, the same idea applies to any holonomic sequence, by deducing the first order matrix recurrence and computing the matrix product using BinSplit.

Now we shall see that the q-holonomic case is similar. First, let q be a positive integer of B bits and consider the computation of the q-factorial u N (q) = (1 + q)(1

+ q + q 2 ) • • • (1 + q + • • • + q N -1 ),
as an illustrative example. For each factor we have the trivial inequalities

q n < 1 + q + • • • + q n < q n+1 ,
meaning that q N (N -1)/2 < u N (q) < q N (N +1)/2 , so the bitsize of u N (q) is of magnitude N 2 B. The "naive" algorithm of deducing u N (q) by first computing the integers q i , then the corresponding sums and products, has Õ(N 3 B) binary complexity. This method is not (quasi-)optimal with respect to the output size. It is also easy to see that the presented baby-steps/giant-steps based algorithms yield bad bit-complexity as well, despite their good arithmetic cost.

Similarly, if

u N (q) = N -1 n=0 q n 2 ,
then the integer u N (q) is bounded in absolute value from above by N q (N -1) 2 and by q (N -1) 2 from below, so its bitsize is again of magnitude N 2 B. The "naive" algorithm consisting of computing the terms q i one after the other before summing, has again non-optimal bit-complexity Õ(N 3 B).

Can one do better? The answer is "yes" and one can even achieve a complexity which is quasi-linear in the bitsize of the output. Similarly to the holonomic setting, it is sufficient to use the q-holonomic character of u N (q), and to reduce its computation to that of a q-factorial matrix as in §2.2, which can then be handled with BinSplit. To be more precise, given an integer or rational number q of bitsize B and any q-holonomic sequence (u n (q)) n≥0 defined by polynomials c 0 , . . . , c r ∈ Z[x, y] with c r (q, q n ) = 0 for any n ∈ N, we define

M (x) ∈ Q(x) by      -cr-1(q,x) cr(q,x) • • • -c1(q,x) cr(q,x) -c0(q,x) cr(q,x) 1 • • • 0 0 . . . . . . . . . . . . 0 • • • 1 0      .
Then, as observed before, u N can be read off from

M (q N -1 ) • • • M (q)M (1),
which we aim to compute efficiently. Again, instead of using baby-steps/giant-steps, it is a better idea to use binary splitting by applying Algorithm 5 to A = [M (1), . . . , M (q N -1 )].

Note that obviously, any element in A is a matrix with rational entries of bitsize bounded by O(N B). Therefore, the complexity of BinSplit does not exceed Õ(N 2 B) by the same argument as before, now using fast multiplication of rational numbers. These considerations prove Theorem 8. Under the assumptions of Theorem 3, with K = Q, the term u N (q) can be computed in Õ(N 2 B) bit operations, where B is the bitsize of q.

As a corollary, (truncated) solutions of q-difference equations can be evaluated using the same (quasi-linear) bit-complexity. This result should be viewed as the q-analogue of the classical fact that holonomic functions can be evaluated fast using binary splitting, a 1988 result by the Chudnovsky brothers [49, §6], anticipated a decade earlier (without proof) by Schroeppel and Salamin in Item 178 of [START_REF] Beeler | Artificial Intelligence Memo[END_REF].

Applications

Combinatorial q-holonomic sequences

As already mentioned, many q-holonomic sequences arise in combinatorics, for example in connection with the enumeration of lattice polygons, where q-analogues of the Catalan numbers 1 n+1 2n n occur naturally [START_REF] Fürlinger | q-Catalan numbers[END_REF][START_REF] Gessel | A noncommutative generalization and q-analog of the Lagrange inversion formula[END_REF], or in the enumeration of special families of matrices with coefficients in the finite field F q [START_REF] Kirillov | On the number of solutions of the equation X 2 = 0 in triangular matrices over a finite field[END_REF][START_REF] Kirillov | On a remarkable sequence of polynomials[END_REF][START_REF] Yip | Rook placements and Jordan forms of upper-triangular nilpotent matrices[END_REF], where sequences related to the Gaussian coefficients n k q also show up.

A huge subfield of combinatorics is the theory of partitions [START_REF] Andrews | The theory of partitions[END_REF], where q-holonomic sequences occur as early as in the famous Rogers-Ramanujan identities [START_REF] Rogers | Proof of certain identities in combinatory analysis[END_REF][START_REF] Rogers | Second Memoir on the Expansion of certain Infinite Products[END_REF], see also [START_REF] Andrews | The theory of partitions[END_REF]Ch. 7], e.g.,

1 + n≥1 q n 2 (1 -q) • • • (1 -q n ) = n≥0 1 (1 -q 5n+1 )(1 -q 5n+4 )
which translates the fact that the number of partitions of n into parts that differ by at least 2 is equal to the number of partitions of n into parts congruent to 1 or 4 modulo 5. Andrews [START_REF] Andrews | Partition identities[END_REF][START_REF] Andrews | A general theory of identities of the Rogers-Ramanujan type[END_REF], see also [START_REF] Andrews | The theory of partitions[END_REF]Chapter 8], laid the foundations of a theory able to capture the q-holonomy of any generating function of a so-called linked partition ideal.

As a consequence, a virtually infinite number of special families of polynomials coming from partitions can be evaluated fast. For instance, the family of truncated polynomials

F n (x) := ∞ k=1 (1 -x k ) 3 mod x n ,
can be evaluated fast due to our results and to the identity [113, §6]

F N (q) = ( n+1 2 )<N
(-1) n (2n + 1)q ( n+1 2 ) .

Evaluation of q-orthogonal polynomials

In the theory of special functions, orthogonal polynomials play a fundamental role. There exists an extension to the q-framework of the theory, see e.g., Chapter 9 in Ernst's book [START_REF] Ernst | A comprehensive treatment of q-calculus[END_REF]. Amongst the most basic examples, the discrete q-Hermite polynomials [START_REF] Al-Salam | Some orthogonal q-polynomials[END_REF][START_REF] Askey | Continuous q-Hermite polynomials when q > 1[END_REF] are defined by their q-exponential generating function n≥0 F n,q (x)

t n [n] q ! = e q (xt) e q (t)e q (-t) , and therefore they satisfy the second-order linear q-recurrence

F n+1,q (x) = xF n,q (x) -(1 -q n )q n-1 F n-1,q (x), n ≥ 1,
with initial conditions F 0,q (x) = 1, F 1,q (x) = x. From there, it follows that for any α ∈ K, the sequence (F n,q (α)) n≥0 is q-holonomic, thus the evaluation of the N -th polynomial at x = α can be computed fast. The same is true for the continuous q-Hermite polynomials, for which 2αH n,q (α) = H n+1,q (α) + (1 -q n )H n-1,q (α) for n ≥ 1, and H 0,q (α) = 1, H 1,q (α) = 2α. More generally, our results in §3 imply that any family of q-orthogonal polynomials can be evaluated fast.

Polynomial and rational solutions of q-difference equations

The computation of polynomial and rational solutions of linear differential equations lies at the heart of several important algorithms, for computing hypergeometric, d'Alembertian and Liouvillian solutions, for factoring and for computing differential Galois groups [START_REF] Abramov | q-hypergeometric solutions of q-difference equations[END_REF][START_REF] Abramov | D'Alembertian solutions of inhomogeneous linear equations (differential, difference, and some other)[END_REF][START_REF] Van Der Put | Galois theory of linear differential equations[END_REF]. Creative telescoping algorithms (of second generation) for multiple 20 integration with parameters [START_REF] Chyzak | An extension of Zeilberger's fast algorithm to general holonomic functions[END_REF][START_REF] Koutschan | A fast approach to creative telescoping[END_REF] also rely on computing rational solutions, or deciding their existence. The situation is completely similar for q-difference equations, i.e. equations of the form Ly = a ν (x)y(q ν x) + a ν-1 (x)y(q ν-1 x)

+ • • • + a 0 (x)y(x) = 0, (11)
with a j (x) ∈ Q[x], for all j = 0, . . . , ν, such that a 0 (x)a ν (x) is not identically zero.

Improving algorithms for polynomial and rational solutions of such equations is important for finding q-hypergeometric solutions [START_REF] Abramov | q-hypergeometric solutions of q-difference equations[END_REF], for computing q-difference Galois groups [START_REF] Arreche | Computing differential Galois groups of second-order linear q-difference equations[END_REF][START_REF] Hendriks | An algorithm for computing a standard form for second-order linear q-difference equations[END_REF], and for performing q-creative telescoping [START_REF] Chyzak | An extension of Zeilberger's fast algorithm to general holonomic functions[END_REF][START_REF] Koornwinder | On Zeilberger's algorithm and its q-analogue[END_REF][START_REF] Koutschan | A fast approach to creative telescoping[END_REF].

In both differential and q-differential cases, algorithms for computing polynomial solutions proceed in two distinct phases: (i) compute a degree bound N , potentially exponentially large in the equation size; (ii) reduce the problem of computing polynomial solutions of degree at most N to linear algebra. Abramov, Bronstein and Petkovšek showed in [START_REF] Abramov | On polynomial solutions of linear operator equations[END_REF] that, in step (ii), linear algebra in size N can be replaced by solving a much smaller system, of polynomial size. However, setting up this smaller system still requires linear time in N , essentially by unrolling a (q-)linear recurrence up to terms of indices close to N . For differential (and difference) equations, this step has been improved in [START_REF] Bostan | Low complexity algorithms for linear recurrences[END_REF][START_REF] Bostan | Fast algorithms for polynomial solutions of linear differential equations[END_REF], by using Chudnovskys' algorithms for computing fast the N -th term of a holonomic sequence. This allows for instance to decide (non-)existence of polynomial solutions in sublinear time Õ( √ N ). Moreover, when polynomial solutions exist, one can represent/manipulate them in compact form using the recurrence and initial terms as a compact data structure. Similar ideas allow to also compute rational solutions in compact form in the same complexity, see [START_REF] Bostan | Algorithmes Efficaces en Calcul Formel[END_REF]Chap. 17].

The same improvements can be transferred to linear q-difference equations, in order to improve the existing algorithms [START_REF] Abramov | On polynomial solutions of linear operator equations[END_REF][START_REF] Abramov | Rational solutions of linear difference and q-difference equations with polynomial coefficients[END_REF][START_REF] Khmel | Improved algorithms for solving difference and q-difference equations[END_REF]. In this case, setting up the smaller system in phase (ii) amounts to computing the N -th term of a q-holonomic sequence, and this can be done fast using our results in §3. A technical subtlety is that, as pointed out in [1, §4.3], it is not obvious in the q-difference case how to guarantee the non-singularity of the q-recurrence on the coefficients of the solution. This induces potential technical complications similar to the ones for polynomial solutions of differential equations in small characteristic, which can nevertheless be overcome by adapting the approach described in [43, §3.2]. Similar improvements can be also transferred to systems [START_REF] Abramov | A direct algorithm to compute rational solutions of first order linear q-difference systems[END_REF][START_REF] Barkatou | Simple forms and rational solutions of pseudo-linear systems[END_REF].

Let us finish this discussion by pointing out briefly an application of these improvements. Desingularizing a linear differential operator L(x, ∂ x ) consists in computing a left multiple with all apparent singularities removed. It is a central task for determining the Weyl closure of L [START_REF] Tsai | Weyl closure of a linear differential operator[END_REF]. The computation of polynomial solutions of Fourier dual operators is a basic step for performing desingularizations [START_REF] Chyzak | Taming apparent singularities via Ore closure[END_REF]. By duality, the order N of the desingularization corresponds to the degree of polynomial solutions of the dual L ⋆ of L. This remark in conjunction with the fast algorithms for polynomial solutions [START_REF] Bostan | Fast algorithms for polynomial solutions of linear differential equations[END_REF], themselves based on the fast computation of matrix factorials, allows to speed up the computation of desingularizations. The situation is similar in other Ore algebras, and in particular for q-difference equations. Therefore, our algorithmic improvements in the computation of polynomial solutions of q-difference equations, themselves based on the fast computation of matrix q-factorials, have a direct impact on the acceleration of the desingularizations process for q-difference equations. It is less obvious to us whether other desingularization algorithms, such as the one from [START_REF] Koutschan | Desingularization in the q-Weyl algebra[END_REF], could also benefit from these remarks.

Computing curvatures of q-difference equations

A natural application of the fast computation of matrix q-factorials is the computation of curvatures of q-difference equations, since in this area these objects appear quite inherently. Another strong motivation comes from the fact that the q-analogue [START_REF] Bézivin | Les suites q-récurrentes linéaires[END_REF] of Grothendieck's conjecture (relating equations over Q with their reductions modulo primes p) is proved [START_REF] Vizio | Arithmetic theory of q-difference equations: the q-analogue of Grothendieck-Katz's conjecture on p-curvatures[END_REF][START_REF] Vizio | Intrinsic approach to Galois theory of q-difference equations, with the preface to Part IV by Anne Granier[END_REF], while the classical differential case is widely open [START_REF] Katz | Algebraic solutions of differential equations (p-curvature and the Hodge filtration)[END_REF]. Still, in the latter setting algorithms have been developed allowing to compute p-curvatures fast [START_REF] Bostan | A fast algorithm for computing the characteristic polynomial of the p-curvature[END_REF][START_REF] Bostan | A fast algorithm for computing the p-curvature[END_REF][START_REF] Bostan | Computation of the similarity class of the p-curvature[END_REF][START_REF] Bostan | Fast algorithms for differential equations in positive characteristic[END_REF]. They allow, for example, to perform a quick heuristic, but reliable in practice, test for the existence of a basis of algebraic solutions of a linear differential operator [START_REF] Bostan | Automatic classification of restricted lattice walks[END_REF].

The q-analogue of Grothendieck's conjecture investigates rational solutions of qdifference equations. Similarly to the differential setting, one naturally associates to the equation ( 11) the linear q-difference system Y (qx) = A(x)Y (x), where

A(x) =      0 1 • • • 0 . . . . . . 0 0 • • • 1 -a 0 /a ν -a 1 /a ν • • • -a ν-1 /a ν      .
Then it easily follows that if z(x) solves [START_REF] Alman | A refined laser method and faster matrix multiplication[END_REF], then (z(x), z(qx), . . . , , z(q ν-1 x)) t is a solution of Y (qx) = A(x)Y (x). Moreover, these equations are in some sense equivalent through the q-analogues of the Wronskian Lemma and the Cyclic Vector Lemma whenever q is not a root of unity or order smaller than ν; see [START_REF] Vizio | Intrinsic approach to Galois theory of q-difference equations, with the preface to Part IV by Anne Granier[END_REF] for details. If q is considered as a variable, then Di Vizio and Hardouin proved that (11) admits a full set of solutions in Q(q, x) if and only if for almost all natural numbers n,

C n (x) := A(q n-1 x) • • • A(qx)A(x) ≡ Id ν mod GL ν (R n (x)),
where R n = Q[q]/Φ n (q), with Φ n (x) the n-th cyclotomic polynomial. The elements in the sequence (C n (x)) n≥1 are known as curvatures of the q-difference system and are clearly just matrix q-factorials.

On the other hand, if q ∈ Q then it already follows from main result of [START_REF] Vizio | Arithmetic theory of q-difference equations: the q-analogue of Grothendieck-Katz's conjecture on p-curvatures[END_REF] that (11) admits a basis of rational solutions in Q(x) if and only if for almost all primes p, A(q κp-1 x) • • • A(qx)A(x) ≡ Id ν mod p ℓ , where κ p = ord p (q) and ℓ p ∈ Z such that 1 -q κp = p ℓp h g , with h, g ∈ Z coprime to p. On these types of questions there is more progress in the q-difference setting than in the classical differential one. Yet, unfortunately, the theorems above are not proven to be effective in the sense that still infinitely many conditions need to be checked in order to conclude the implication we are mostly interested in. Therefore, the computation of any finite number of curvatures only provides a heuristic for the existence of rational solutions of a q-difference equation. Moreover, the mentioned algorithms in §4.3 compute rational solutions of equations of type [START_REF] Alman | A refined laser method and faster matrix multiplication[END_REF] and therefore allow to decide rigorously about the existence of such a basis. However, all these methods have a cost which is potentially exponential in the size of the input. Our goal in this section is to design a fast heuristic test for the existence of a basis of rational solutions of a q-difference equation using curvatures.

If q is a variable, we want to check whether C n (x) ≡ Id ν mod GL ν (R n (x)) for many n. Clearly, after the reduction mod Φ n (q), the polynomial C n (x) has arithmetic size n 2 over Q and Õ(n 3 ) bitsize. Hence, computing C n (x) is unnecessarily costly. We propose to work over R n,p := F p [q]/Φ n (q) for some (large) prime p; moreover, we compute C n (x 0 ) for some randomly chosen x 0 ∈ F p . Furthermore, in order to avoid computing general cyclotomic polynomials, we compute C n (x) only for prime numbers n. After these considerations, it is easy to see that Algorithm 3 applies and allows to deduce

C n (x 0 ) modulo GL ν (R n,p ) in Õ( √ n) arithmetic operations in R n,p , hence Õ(n 3/2 ) operations in F p .
Finally, if we want to deduce this quantity for all primes n between 2 and some N ∈ N, it is wiser to apply the accumulating remainder tree method presented in [START_REF] Costa | A search for Wilson primes[END_REF][START_REF] Harvey | Counting points on hyperelliptic curves in average polynomial time[END_REF], which allows for quasi-optimal complexity of Õ(N 2 ) in this case.

If q is some rational number and the goal is to test whether Y (qx) = A(x)Y (x) has a full set of rational solutions, one may check

C κp (x) = A(q κp-1 x) • • • A(qx)A(x) ≡ Id ν mod p ℓ ,
for many primes p. Unfortunately, finding the order κ p in practice may be costly, therefore we shall check the weaker assumption C p-1 (x) ≡ Id ν mod p. Conjecturally, this equality for sufficiently many primes p is also enough to conclude on the existence of a basis of rational solutions. The presented Algorithm 3 allows to compute C p-1 (x 0 ) mod p for some x 0 ∈ F p in Õ( √ p) arithmetic cost. Finally, again, if we choose N, x 0 ∈ N, then the accumulating remainder tree allows to deduce C p-1 (x 0 ) mod p for all primes p between 2 and N quasi-optimally in Õ(N ) bit operations.

q-hypergeometric creative telescoping

In the case of differential and difference hypergeometric creative telescoping, it was demonstrated in [START_REF] Bostan | Low complexity algorithms for linear recurrences[END_REF] that the compact representation for polynomial solutions can be used as an efficient data structure, and can be applied to speed up the computation of Gosper forms and Zeilberger's classical summation algorithm [119, §6]. The key to these improvements lies in the fast computation of the N -th term of a holonomic sequence, together with the close relation between Gosper's algorithm and the algorithms for rational solutions.

Similarly, in the q-difference case, Koornwinder's q-Gosper algorithm [97, §5] is closely connected to Abramov's algorithm for computing rational solutions [2, §2], and this makes it possible to transfer the improvements for rational solutions to the q-Gosper algorithm. This leads in turn to improvements upon Koornwinder's algorithm for q-hypergeometric summation [START_REF] Koornwinder | On Zeilberger's algorithm and its q-analogue[END_REF], along the same lines as in the differential and difference cases [START_REF] Bostan | Low complexity algorithms for linear recurrences[END_REF].

Experiments

Algorithms 1 and 2 were implemented in Magma and Algorithm 3 in Maple. All implementations deliver some encouraging timings. Of course, since these algorithms are designed to be fast in the arithmetic model, it is natural to make experiments over a finite field K, or over truncations of real/complex numbers, as was done in [START_REF] Nogneng | On the evaluation of some sparse polynomials[END_REF] for the problem in §2.2.

Recall that both Algorithms 1 and 2 compute N -1 i=0 (α -q i ) ∈ K, given α, q in a field K, and N ∈ N, whereas Algorithm 3 finds M (q N -1 ) • • • M (q)M (1) ∈ K n×n for a degree N Naive algorithm Algorithm 1 Algorithm 2

3355.39

Table 1 Comparative timings (in seconds) for the computation of N-1 i=0 (αq i ) ∈ Fp, with p = 2 30 + 3 and (α, q) randomly chosen in Fp × Fp. All algorithms were executed on the same machine, running Magma v. 2.24. For each target degree N , each execution was limited to 1 hour. Naive algorithm could reach degree N = 2 32 , Algorithm 1 degree N = 2 46 given polynomial matrix M (x) ∈ K[x] of size n×n and q ∈ K, N ∈ N. In our experiments, K is the finite field F p with p = 2 30 + 3 elements.

Timings for Algorithms 1 and 2 are presented in Table 1. We compare the straightforward iterative algorithm (column Naive), to the fast baby-step/giant-step algorithms, one based on subproduct trees and resultants (column Algorithm 1), the other based on multipoint evaluation on geometric sequences (column Algorithm 2).

Some conclusions can be drawn by analyzing these timings:

• The theoretical complexities are perfectly reflected in practice: as N is increased from 2 2k to 2 2k+2 , timings are also multiplied (roughly) by 4 in column Naive, and (roughly) by 2 in columns Algorithm 1 and Algorithm 2.

• The asymptotic regime is reached from the very beginning.

• Algorithm 2 is always faster than Algorithm 1, which is itself much faster than the Naive algorithm, as expected.

• A closer look into the timings shows that for Algorithm 1, ≈ 80% of the time is spent in step (3) (resultant computation), the other steps taking ≈ 10% each; for Algorithm 2, step (1) takes ≈ 25%, step (2) takes ≈ 75%, and step (3) is negligible. 24 In order to visualize the performance of Algorithm 3, we took random polynomial matrices in M n (F p [x]) of degree d. Figure 1 compares the time needed to compute the Matrix q-factorial for d = 1 and n = 2, 4, 8 as N growths. The black lines represent the best linear fits to the data points and are given by the equations y = 0.56x -14.1, y = 0.59x -12.5 and y = 0.53x -9.0 respectively. Note that we are plotting on a log-log scale, therefore the established complexity of Õ(N 1/2 ) is indicated by the coefficient of x in these linear fits which is always only slightly greater than 1/2. In the same figure, we also show the timings for d = 1 and n = 2 of the naive algorithm given by successively computing and multiplying M (q i ) together. The best linear fit almost perfectly describes this data and has a slope of 1.04 ≈ 1, in line with the linear complexity in N .

2
In Figure 2 we show similar timings, however now for n = 2 and d = 2, 4, 8. The linear fits to the data are now given by y = 0.54x -12.7, y = 0.54x -11.9 and y = 0.56x -11.5. Again, they describe the observations very well and the coefficients of the regressions are in line with the proven complexity. We observe that, as expected, the lines have slopes of roughly 1/2 and are closer together than in the previous figure. The naive method has a slope of 1.02 ≈ 1.

Conclusion and future work

We have shown that selected terms of q-holonomic sequences can be computed fast, both in theory and in practice, the key being the extension of classical algorithms in the holonomic ("q = 1") case. We have demonstrated through several examples that this basic algorithmic improvement has many other algorithmic implications, notably on the faster evaluation of many families of polynomials and on the acceleration of algorithms for q-difference equations.

Here are some questions that should be investigated in the future. modulo p, for one [START_REF] Bostan | Linear recurrences with polynomial coefficients and application to integer factorization and Cartier-Manin operator[END_REF] or several [START_REF] Harvey | Counting points on hyperelliptic curves in average polynomial time[END_REF] primes p. A natural extension is to ask the same with G p (x) replaced by p-1 2 k=1 (q 3k x 3 + aq k x + b). This might have applications related to §4.4, or to counting points on q-deformations [START_REF] Scholze | Canonical q-deformations in arithmetic geometry[END_REF].

2. (Computing q-deformed real numbers) Recently, Morier-Genoud and Ovsienko [START_REF] Morier-Genoud | On q-deformed real numbers[END_REF] introduced q-analogues of real numbers, see also [START_REF] Stum | On quantum integers and rationals[END_REF][START_REF] Morier-Genoud | q-deformed rationals and q-continued fractions[END_REF]. How fast can one compute (truncations or evaluations of) quantized versions of numbers like e or π?

3. (Evaluating more polynomials) Is it possible to evaluate fast polynomials of the form N ℓ=0 x ℓ s , for s ≥ 3, and many others that escape the q-holonomic class? E.g., [START_REF] Bernstein | Faster computation of isogenies of large prime degree[END_REF] presents a beautiful generalization of Algorithm 1 to the fast evaluation of isogenies between elliptic curves, by using elliptic resultants, with applications to isogeny-based cryptography. 4. ("Precise" complexity of q-holonomic sequences) Can one prove non-trivial lower bounds, ideally matching the upper bounds, on examples treated in this article?

  (p,a,b) N (q), the dependence in a, b of Algorithm 3 is quite high (quasi-linear in a and b). If a and b are fixed and considered as O(1) this dependence is invisible, but otherwise the following variant has the same complexity with respect to N , and a much better cost with respect to a and b. It is based on the simple observation that, if M (x) denotes the polynomial matrix

  indices; moreover d 1 := ⌈d 0 /2⌉. We iterate this process at most ℓ := ⌈log(d 0 )⌉ many times until d ℓ = 1. Then we can easily find U N1 , . . . , U Nn , from which we finally deduce u N1 , . . . , u Nn .Each such step has a cost of at most O(M(n)) = O(M( √ N )) base field operations. Moreover, after first employing Algorithm 3 and by the consideration above, we compute U 1 , U s , . . . , U s 2 in O(M( √ N )) base operations (Theorem 5). Hence, we may choose d 0 = s and for each j = 1, . . . , n let k (0) j be the largest element in {1, s, . . . , (s -1)s} such that k (0) j ≤ N j . Clearly, all conditions of the above refinement step are satisfied and we need at most ⌈log(s)⌉ = O(log N ) many such steps. The total complexity is henceforth

  , and Algorithm 2 degree N = 2 54 = 8 014 398 509 481 984. By extrapolation, the Naive algorithm would have needed ≈ 4 11 × 3028.25 sec. ≈ 400 years on the same instance, and Algorithm 2 approximately 18 hours.

2 8 N

 28 Time in secondsd = 1; n = 2; naive d = 1; n = 2 d = 1; n = 4 d = 1; n = 8

Figure 1 :

 1 Figure 1: Timings of Algorithm 3 implemented in Maple 2020.2. We compare d = 1 and n = 2, 4, 8 for various values of N .

2 10 2 8 N

 28 14 2 18 2 22 2 26 2 30 2 34 2 38 Time in seconds d = n = 2; naive d = 2; n = 2 d = 4; n = 2 d = 8; n = 2

Figure 2 :

 2 Figure 2: Timings of Algorithm 3 implemented in Maple 2020.2. We compare n = 2 and d = 2, 4, 8 for various values of N .

  polynomial multiplication

	10:	else
	11:

p 1 (x) ← BS((t -1)/2) 12:

  10 2 14 2 18 2 22 2 26 2 30 2 34 2 38

	2 6
	2 4
	2 2
	2 0
	2 -2
	2 -4
	2 -6

If q n = 1 for some n < N , then it is enough to compute the product of α -q i for i = 0, . . . , n -1 and its appropriate power. The latter step can be done efficiently (in essentially log(N ) operations) using binary powering. Our main interest lies therefore in q ∈ K that are not roots of unity of small order.

Acknowledgements. We would like to thank Luca De Feo for his initial question, who motivated this work, and for the very interesting subsequent discussions. Our warm thanks go to Lucia Di Vizio for valuable comments that led to §4.4. Moreover, we thank her and Kilian Raschel for their careful reading of the first version of this manuscript. The first author was supported in part by DeRerumNatura ANR-19-CE40-0018 and the second author by the Austrian Science Fund (FWF) P-31338.