Fast Computation of the N -th Term of a q -Holonomic Sequence and Applications
aa r X i v : . [ c s . S C ] D ec Fast Computation of the N -th Term of a q -Holonomic Sequenceand Applications Alin Bostan a , Sergey Yurkevich b a Inria, Université Paris-Saclay (France) b University of Vienna (Austria)
Abstract
In 1977, Strassen invented a famous baby-step/giant-step algorithm that computes thefactorial N ! in arithmetic complexity quasi-linear in √ N . In 1988, the Chudnovskybrothers generalized Strassen’s algorithm to the computation of the N -th term of anyholonomic sequence in essentially the same arithmetic complexity. We design q -analoguesof these algorithms. We first extend Strassen’s algorithm to the computation of the q -factorial of N , then Chudnovskys’ algorithm to the computation of the N -th term ofany q -holonomic sequence. Both algorithms work in arithmetic complexity quasi-linearin √ N ; surprisingly, they are simpler than their analogues in the holonomic case. Weprovide a detailed cost analysis, in both arithmetic and bit complexity models. Moreover,we describe various algorithmic consequences, including the acceleration of polynomialand rational solving of linear q -differential equations, and the fast evaluation of largeclasses of polynomials, including a family recently considered by Nogneng and Schost.
1. Introduction
A classical question in algebraic complexity theory is: how fast can one evaluate aunivariate polynomial at one point? The precise formulation of this question depends onthe model of computation. We will mainly focus on the arithmetic complexity model, inwhich one counts base field operations at unit cost.Horner’s rule evaluates a polynomial P in O (deg( P )) operations. Ostrowski [112]conjectured in 1954 that this is optimal for generic polynomials , i.e. , whose coefficientsare algebraically independent over the prime subfield. This optimality result was proveda few years later by Pan [114].However, most polynomials that one might wish to evaluate “in practice” have coeffi-cients which are not algebraically independent. Paterson and Stockmeyer [116] 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 algorithmuses a linear amount of scalar multiplications, so it is not well adapted to the evaluation Email addresses: [email protected] (Alin Bostan), [email protected] (SergeyYurkevich)
Preprint submitted to Elsevier December 17, 2020 t points from the base field K , since in this case the total arithmetic complexity, countedin 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 − + · · · + 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, afamily ( 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 ) , whateverthe algorithm used.Paterson and Stockmeyer [115, 116] proved the existence of polynomials in Q [ x ] whichare hard to compute (note that this does not follow from Pan’s result [114]). However,their proof was based on a non-constructive argument. Specific families of hard-to-compute polynomials were first exhibited by Strassen [133]. For instance, he provedthat for large N , the polynomial P Nℓ =0 ℓ x ℓ needs at least p N/ (3 log N ) operationsto be evaluated. The techniques were refined and improved by Borodin and Cook [29],Lipton [104] and Schnorr [127], who produced explicit examples of degree- N polynomialswhose evaluation requires a number of operations linear in √ N . Subsequently, variousmethods have been developed to produce similar results on lower bounds , e.g. , by Heintzand Sieveking [83] using algebraic geometry, and by Aldaz et al. [10] using a combinatorialapproach. The topic is vast and very well summarized in the book by Bürgisser, Clausenand Shokrollahi [45].In this article, we focus on upper bounds , that is on the design of fast algorithms forspecial families of polynomials, which are hard to compute, but easier to evaluate thangeneric polynomials. For instance, for the degree- (cid:0) N (cid:1) polynomial Q N ( x ) := P ( 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 ). Anotherstriking example is R N ( x ) := N X ℓ =0 x ℓ , of degree N , and whose evaluation can also be performed in complexity quasi-linearin √ N , as shown recently by Nogneng and Schost [110] (see §2.2). In both cases, thesecomplexities are obtained by clever although somehow ad-hoc algorithms. The startingpoint of our work was the question whether these algorithms for Q N ( x ) and R N ( x ) couldbe treated in a unified way, which would allow to evaluate other families of polynomialsin a similar complexity.The answer to this question turns out to be positive. The key idea, very simple andnatural, 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 recur-rence with polynomial coefficients in q and q n , how fast can one compute its N -th term? 2n the more classical case of holonomic sequences (satisfying linear recurrences withpolynomial 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 [134].The Chudnovsky brothers extended in [49] Strassen’s algorithm to the computation ofthe 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 sim-pler in the q -holonomic case than in the usual holonomic setting, essentially becausemultipoint evaluation on arithmetic progressions used as a subroutine in Strassen’s andChudnovskys’ algorithms is replaced by multipoint evaluation on geometric progressions,which is considerably simpler [42].A consequence of our results is that the following apparently unrelated polynomialsand rational functions can be evaluated fast (note the change in notation, with thevariable x denoted now by q ):• A n ( q ) , the generating function of the number of partitions into n positive integerseach occurring at most twice [141], i.e. , the coefficient of t n in the product Y k ≥ (1 + q k t + q k t ) . • B n ( q ) := Q ∞ i =1 (1 − q i ) mod q n ; by Euler’s pentagonal theorem [113, §5], B n ( q ) = 1 + X i (3 i +1) < n ( − i (cid:16) q i (3 i − + q i (3 i +1)2 (cid:17) . • The number C n ( q ) of n × n upper-triangular matrices over F q (the finite fieldwith q elements), whose square is the zero matrix [94]; by [59], C n ( q ) is equal to C n ( q ) = X j (cid:20)(cid:18) nn − j (cid:19) − (cid:18) nn − j − (cid:19)(cid:21) · q n − j − j . The common feature, exploited by the new algorithm, is that the sequences ( A n ( q )) n ≥ , ( B n ( q )) n ≥ , ( C n ( q )) n ≥ are all q -holonomic. Actually, q -holonomic sequences are ubiq-uitous, so the range of application of our results is quite broad. This stems from thefact that they are coefficient sequences of power series satisfying q -difference equations,or equivalently, q -shift (or, q -differential) equations. From that perspective, our topicbecomes intimately connected with q -calculus. The roots of q -calculus are in works offamous mathematicians such as Rothe [125], Gauss [75] and Heine [82]. The topic gainedrenewed interest in the first half of the 20th century, with the work, both on the formaland analytic aspects, of Tanner [135], Jackson [87–89], Carmichael [47], Mason [106],Adams [7, 8], Trjitzinsky [137], Le Caine [102] and Hahn [77], to name just a few. Mod-ern accounts of the various aspects of the theory (including historical ones) can be foundin [57, 60, 95].One of the reasons for interest in q -differential equations is that, formally, as q tendsto , the q -derivative f ( qx ) − f ( x )( q − x tends to f ′ ( x ) , thus to every differential equation corre-sponds a q -differential equation which goes formally to the differential equation as q → .3n nice cases, (some of) the solutions of the q -difference equation go to solutions of theassociated differential equation as q → . An early example of such a good deformationbehavior is given by the basic hypergeometric equation of Heine [82], see also [95, §1.10].In computer algebra, q -holonomic sequences were considered starting from the earlynineties, in the context of computer-generated proofs of identities in the seminal paperby Wilf and Zeilberger [140], notably in Section 5 (“Generalization to q -sums and q -multisums”) and in Section 6.4 (“ q -sums and integrals”). Creative telescoping algorithmsfor (proper) q -hypergeometric sequences are discussed in various references [28, 48, 119];several implementations of those algorithms are described for instance in [91, 118, 122,132]. Algorithms for computing polynomial, rational and q -hypergeometric solutions of q -difference equations were designed by Abramov and collaborators [1, 2, 4, 92]. Thesealgorithms are important for several reasons. One is that they lie at the heart of the vastgeneralization by Chyzak [50, 51] of the Wilf and Zeilberger algorithmic theory, for thetreatment of general q -holonomic (not only q -hypergeometric) symbolic summation andintegration via creative telescoping. In that context, a multivariate notion of q -holonomyis needed; the foundations of the theory were laid by Zeilberger [144] and Sabbah [126](in the language of D-modules), see also [48, § 2.5] and [72].The simplest non-trivial holonomic sequence is ( n !) n ≥ , whose n -th term combinato-rially counts the number of permutations of n objects. If instead of direct counting, oneassigns to every permutation π its number of inversions inv( π ) , i.e. , the number of pairs ≤ i < j ≤ n with π ( i ) > π ( j ) , the refined count (by size and number of inversions) is [ n ] q ! := (1 + q )(1 + q + q ) · · · (1 + q + · · · + q n − ) . 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 Gaussiancoefficients , defined by (cid:18) nk (cid:19) q := [ n ] q ![ k ] q ![ n − k ] q ! . They have many counting interpretations, e.g., they count the k -dimensional subspacesof F nq (points on Grassmannians over F q ). There are q -analogs to (almost) everything. Toselect just two more basic examples, the q -analog [14, Thm. 3.3] of the binomial theoremis given by n Y k =1 (1 + q k − x ) = n X k =0 (cid:18) nk (cid:19) q q ( k ) x k , (1)and the q -version [14, Thm. 3.4] of the Chu-Vandermonde identity is n X k =0 q k (cid:18) mk (cid:19) q (cid:18) nk (cid:19) q = (cid:18) m + nn (cid:19) q . (2)The ubiquity of q -holonomic sequences is manifest in plenty of fields: partitiontheory [14, 15, 105, 113, 131, 141] and other subfields of combinatorics [16, 44, 59,63, 93, 94, 142]; theta functions and modular forms [22, 73, 100, 101, 143]; specialfunctions [30, 86, 95] and in particular orthogonal polynomials [97]; algebraic geom-etry [58], representation theory [85]; knot theory [68–72]; Galois theory [84]; numbertheory [6, 55, 111]. 4he main messages of this article are that for any example of a q -holonomicsequence occurring in those various fields, one can compute selected coeffi-cients faster than by a direct algorithm and that this fact finds a tremendousnumber of applications . Complexity basics.
We estimate the arithmetic complexities of algorithms by count-ing arithmetic operations (+ , − , × , ÷ ) in the base field K at unit cost. We use stan-dard 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 upperbound is θ < . [11, 64]. As usual, O ( · ) stands for the big-Oh notation and ˜ O ( · ) is used to hide polylogarithmic factors in the argument. Most arithmetic operations onunivariate polynomials of degree d in K [ x ] can be performed in quasi-linear complex-ity ˜ O ( d ) : multiplication, shift, interpolation, gcd, resultant, etc . A key feature of theseresults is the reduction to fast polynomial multiplication, which can be performed intime M ( d ) = O ( d log d log log d ) [46, 129]. Finally, the arithmetic cost of multiplicationof polynomial matrices of size n and degree d is denoted by MM ( n, d ) and we have MM ( n, d ) = O ( n θ d + n M ( d )) = ˜ O ( n θ d ) [42]. An excellent general reference for thesequestions is the book by von zur Gathen and Gerhard [74].A short version of this article has appeared at the ISSAC’20 conference [31]. In thepresent version, we included the proofs of Theorems 6 and 8, we added a new Theo-rem 5 containing a detailed complexity analysis of the main algorithm ( Algorithm 3 ) withrespect to all parameters, and we displayed pseudo-code for the algorithms as well asfigures visualizing their performance. We also elaborated on a task which was mentionedas future work in the previous version, namely the application of our methods to thecomputation 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 evalu-ating Q N ( x ) and R N ( x ) . We show that these are two instances of the same problem andprovide Algorithm 3 which solves both in O ( M ( √ N )) arithmetic complexity. Section 3 isdevoted to the main results; we prove there that Algorithm 3 can be used for computingterms of any q -holonomic sequence with the same cost, and provide extensions and moreinsight. In the same section we also consider the bit-complexity model. We identify andelaborate on several applications for our result in Section 4. In Section 5 we report onimplementations of our algorithms, which deliver encouraging timings, and we finallydescribe future tasks and investigation fields in Section 6.
2. Two motivating examples
Before presenting our main results in Section 3, we describe in this section the ap-proach and main ideas on two basic examples. Both examples concern the fast evaluationof special families of univariate polynomials. In §2.1, we consider polynomials of the form Q ℓ ( x − q ℓ ) , and in §2.2 sparse polynomials of the form P ℓ p ℓ x aℓ + bℓ . In both cases, wefirst 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.5 .1. Evaluation of some structured polynomials Here is our first example, that emerged from a question asked to the first author byLuca De Feo (private email communication, 10 January 2020); this was the starting pointof the article.Let q be an element of the field K , and consider the polynomial F ( x ) := N − Y 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 binarypowering. We assume in what follows that q is nonzero. Obviously, a direct algorithmconsists in computing the successive powers q, q , . . . , q N − using O ( N ) operations in K ,then computing the elements α − , α − q, . . . , α − q N − in O ( N ) more operations in K ,and finally returning their product. The total arithmetic cost of this algorithm 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 baby-step/giant-step strategy, in which, in order to simplify things, we assume that N is aperfect square, N = s . Algorithm 1
1. (Baby-step) Compute the values of q, q , . . . , q s − , and deduce the coefficients ofthe polynomial G ( x ) := s − Y j =0 ( x − q j ) .
2. (Giant-step) Compute Q := q s , Q , . . . , Q s − , and deduce the coefficients of thepolynomial H ( x ) := s − Y k =0 ( α − Q k · x ) .
3. Return the resultant Res ( G, H ) .By the basic property of resultants, the output of this algorithm isRes ( G, H ) = s − Y j =0 H ( q j ) = s − Y j =0 s − Y k =0 (cid:0) α − q sk + j (cid:1) = N − Y i =0 ( α − q i ) = F ( α ) . Using the fast subproduct tree algorithm [74, Algorithm 10.3], one can perform thebaby-step (1) as well as the giant-step (2) in O ( M ( √ N ) log N ) operations in K , and If q n = 1 for some n < N , then it is enough to compute the product of α − q i for i = 0 , . . . , n − andits appropriate power. The latter step can be done efficiently (in essentially log( N ) operations) usingbinary powering. Our main interest lies therefore in q ∈ K that are not roots of unity of small order.
6y [74, Corollary 11.19] the same cost can be achieved for the resultant computation instep (3). Using fast polynomial multiplication, we conclude that F ( α ) can be computedin arithmetic complexity quasi-linear in √ N .Note that if N is not a perfect square, then one can compute F ( α ) as F ( α ) = F ( α ) F ( α ) , where F ( α ) := Q ⌊√ N ⌋ − i =0 ( α − q i ) is computed as in Algorithm 1 , while F ( α ) := Q N − i = ⌊√ N ⌋ ( α − q i ) can be computed naively, since N − ⌊√ N ⌋ = O ( √ N ) .It is possible to speed up the previous algorithm by a logarithmic factor in N usinga slightly different scheme, still based on a baby-step/giant-step strategy, but exploitingthe fact that the roots of F are in geometric progression. Again, we assume that N = s is a perfect square. This alternative algorithm goes as follows. Note that it is very closein spirit to Pollard’s algorithm described on page 523 of [120]. Algorithm 2
1. (Baby-step) Compute q, q , . . . , q s − , and deduce the coefficients of the polynomial P ( x ) := Q s − j =0 ( α − q j · x ) .2. (Giant-step) Compute Q := q s , Q , . . . , Q s − , and evaluate P simultaneously at , Q, . . . , Q s − .3. Return the product P ( Q s − ) · · · P ( Q ) P (1) .Obviously, the output of this algorithm is s − Y k =0 P ( Q k ) = s − Y k =0 s − Y j =0 ( α − q j · q sk ) = N − Y 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 ) = Q s − j =0 ( α − q j · x ) in step (1) without computing the subproduct tree, by usinga divide-and-conquer scheme which exploits the fact that P t ( x ) = P t ( q t x ) · P t ( x ) and P t +1 ( x ) = ( α − q t x ) · P t ( q t x ) · P t ( x ) . The cost of this algorithm is O ( M ( √ N )) operationsin K .As for step (2), one can use the fast chirp transform algorithms of Rabiner, Schaferand Rader [121] and of Bluestein [27]. These algorithms rely on the following observation:writing Q ij = Q ( i + j ) · Q − ( i ) · Q − ( j ) and P ( x ) = P sj =0 c j x j implies that the needed values P ( Q i ) = P sj =0 c j Q ij , ≤ i < s , are P ( Q i ) = Q − ( i ) · s X j =0 c j Q − ( j ) · Q ( i + j ) , ≤ i < s, in which the sum is simply the coefficient of x s + i in the product s X j =0 c j Q − ( j ) x s − j s X ℓ =0 Q ( ℓ ) x ℓ ! . This polynomial product can be computed in M ( s ) operations (and even in M ( s )+ O ( s ) using the transposition principle [40, 78], since only the median coefficients x s , . . . , x s − O ( M ( √ N )) oper-ations 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 becauseit is faster by a factor log( N ) compared to the first one, but more importantly because ithas 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 willdo so for the more general case (Algorithm 3). Let us now consider the sequence of sparse polynomial sums v ( p,a,b ) N ( q ) = N − X n =0 p n q an + bn , where p ∈ K and a, b ∈ Q such that a, a + b are both integers. Typical examples are(truncated) modular forms [117], which are ubiquitous in complex analysis [22], numbertheory [143] and combinatorics [14]. For instance, the Jacobi theta function ϑ dependson two complex variables z ∈ C , and τ ∈ C with ℑ ( τ ) > , and it is defined by ϑ ( z ; τ ) = ∞ X n = −∞ e πi ( n τ +2 nz ) = 1 + 2 ∞ X n =1 η n q n , where q = e πiτ is the nome ( | q | < ) and η = e πiz . Here, K = C . Another exampleis the Dedekind eta function , appearing in Euler’s famous pentagonal theorem [113, §5],which has a similar form q · ∞ X n =1 ( − n (cid:16) q n (3 n − + q n (3 n +1)2 (cid:17)! , with q = e πiτ . Moreover, sums of the form v (1 ,a,b ) N ( q ) = P N − n =0 q an + bn , over K = Q or K = F , cruciallyoccur in a recent algorithm by Tao, Crott and Helfgott [136] for the efficient constructionof prime numbers in given intervals, e.g. , in the context of effective versions of Bertrand’spostulate. Actually, (the proof of) Lemma 3.1 in [136] contains the first sublinear com-plexity result for the evaluation of the sum v ( p,a,b ) N ( q ) at an arbitrary point q ; namely,the cost is O ( N θ/ ) , where θ ∈ [2 , is any feasible exponent for matrix multiplication.Subsequently, Nogneng and Schost [110] designed a faster algorithm, and lowered thecost down to ˜ O ( √ 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 [110]. Assume asbefore that N is a perfect square, N = s . The starting point is the remark that v ( p,a,b ) N ( q ) = N − X n =0 p n q an + bn = s − X k =0 s − X j =0 p j + sk q a ( j + sk ) + b ( j + sk ) can be written s − X k =0 p sk q as k + bsk · P ( q ask ) , where P ( y ) := s − X j =0 p j q aj + bj y j . v ( p,a,b ) N ( q ) can be reduced essentially to the simultaneousevaluation of the polynomial P at s = 1 + deg( P ) points (in geometric progression), witharithmetic cost O ( M ( √ N )) .We now describe an alternative algorithm, of similar complexity O ( M ( √ N )) , with aslightly larger constant in the big-Oh estimate, but whose advantage is its potential ofgenerality.Let us denote by u n ( q ) the summand p n q an + bn . Clearly, the sequence ( u n ( q )) n ≥ satisfies the recurrence relation u n +1 ( q ) = A ( q, q n ) · u n ( q ) , where A ( x, y ) := px a + b y a . As an immediate consequence, the sequence with general term v n ( q ) := P n − 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 ( q ) = 0 and v ( q ) = 1 . This scalar recurrence of order two isequivalent to the first-order matrix recurrence (cid:20) v n +2 v n +1 (cid:21) = (cid:20) A ( q, q n ) + 1 − A ( q, q n )1 0 (cid:21) × (cid:20) v n +1 v n (cid:21) . By unrolling this matrix recurrence, we deduce that (cid:20) v n +1 v n (cid:21) = M ( q n − ) (cid:20) v n v n − (cid:21) = M ( q n − ) · · · M ( q ) M (1) × (cid:20) (cid:21) , where M ( x ) := (cid:20) pq a + b x a + 1 − pq a + b x a (cid:21) , hence v N = (cid:2) (cid:3) × M ( q N − ) · · · M ( q ) M (1) × (cid:20) (cid:21) . Therefore, the computation of v N reduces to the computation of the “matrix q -factorial” M ( q N − ) · · · M ( q ) M (1) , whichcan be performed fast by using a baby-step/giant-step strategy similar to the one of thesecond algorithm in §2.1. Again, we assume for simplicity that N = s is a perfectsquare. The algorithm goes as follows. Algorithm 3 (matrix q -factorial)(1) (Baby-step) Compute q, q , . . . , q s − ; deduce the coefficients of the polynomial ma-trix P ( x ) := M ( q s − x ) · · · M ( qx ) M ( x ) .(2) (Giant-step) Compute Q := q s , Q , . . . , Q s − , and evaluate (the entries of) P ( x ) simultaneously at , Q, . . . , Q s − .(3) Return the product P ( Q s − ) · · · P ( Q ) P (1) .Clearly, this algorithm generalizes Algorithm 2 in §2.1 and, as promised, we also providea detailed pseudo-code implementation:
Step1 and
Step2 & Step3 :9 lgorithm 3 ( Step1 ) Input : s, q, M ( x ) Output : M ( q s − x ) · · · M ( qx ) M ( x ) q s ← [ q, q , . . . , q s − ] t ← s function BS ( t ) if t = 1 then return M ( x ) if t is even then p ( x ) ← BS ( t/ p ( x ) ← p ( q t/ x ) ⊲ Using q s return p ( x ) · p ( x ) ⊲ Fast polynomial multiplication else p ( x ) ← BS (( t − / p ( x ) ← p ( q ( t − / x ) ⊲ Using q s p ( x ) ← M ( q t − x ) ⊲ Using q s return p ( x ) · p ( x ) · p ( x ) ⊲ Fast polynomial multiplication
Algorithm 3 ( Step2 & Step3 ) Input : s, Q, P ( x ) Output : P ( Q s − ) · · · P (1) Assumptions: Q = 0 , P ( x ) polynomial matrix of size n × n and degree d ≥ s . Q d ← [ Q, Q , . . . , Q d − ] Q ′ ← /Q Q ′ d ← [ Q − ( d ) , . . . , Q − ( ) , Q ( ) , . . . , Q ( d )] ⊲ Using Q d and Q ′ P s − , . . . , P ⊲ Empty n × n matrices for i from 1 to n do for j from 1 to n do p ( x ) ← P ( x ) i,j ⊲ p ( x ) = c + c x + · · · + c d x d p ( x ) ← P dℓ =0 c ℓ Q − ( ℓ ) x d − ℓ ⊲ Using Q ′ d p ( x ) ← P dℓ =0 Q ( ℓ ) x ℓ ⊲ Using Q ′ d p ( x ) ← P ( x ) · P ( x ) ⊲ p ( x ) = P dℓ =0 r ℓ x ℓ ; fast multiplication for k from 0 to s − do ( P k ) i,j ← r d + k ⊲ P ℓ = P ( Q ℓ ) for ℓ = 0 , . . . , s − P ← for k from to s − do 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 notonly with respect to N , but to the matrix size and degree as well.We remark that when applied to the computation of v ( p,a,b ) N ( q ) , the dependence in a, b Algorithm 3 is quite high (quasi-linear in a and b ). If a and b are fixed and consideredas O (1) this dependence is invisible, but otherwise the following variant has the samecomplexity with respect to N , and a much better cost with respect to a and b . It is basedon the simple observation that, if ˜ M ( x ) denotes the polynomial matrix ˜ M ( x ) := (cid:20) prx + 1 − prx (cid:21) , with r := q a + b , (5)and if ˜ q := q a , then the following matrix q -factorials coincide: M ( q N − ) · · · M ( q ) M (1) = ˜ M (˜ q N − ) · · · ˜ M (˜ q ) ˜ M (1) . Algorithm 4 (matrix q -factorial, variant)(0) (Precomputation) Compute r := q a + b , ˜ q := q a , and ˜ M in (5).(1) (Baby-step) Compute ˜ q, ˜ q , . . . , ˜ q s − ; deduce the coefficients of the polynomial ma-trix ˜ P ( x ) := ˜ M (˜ q s − x ) · · · ˜ M (˜ qx ) ˜ M ( x ) . (2) (Giant-step) Compute ˜ Q := ˜ q s , ˜ Q , . . . , ˜ Q s − , and evaluate (the entries of) ˜ P ( x ) simultaneously at , ˜ Q, . . . , ˜ Q s − .(3) Return the product ˜ P ( ˜ Q s − ) · · · ˜ P ( ˜ Q ) ˜ P (1) .Using binary powering, the cost of the additional precomputation in step (0) is onlylogarithmic in a and b . In exchange, the new steps (2) and (3) are performed on matriceswhose degrees do not depend on a and b anymore (in the previous, unoptimized, versionthe degrees of the polynomial matrices were linear in a and b ). The total arithmetic costwith respect to N is still quasi-linear in √ N .In the next section, we will show that Algorithm 3 can be employed for the fastcomputation 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. Inany case, we can improve on the quasi-linear cost in the degree d of the polynomialmatrix M ( x ) in Algorithm 3 , obtaining a complexity of essentially √ d ; in essence, theidea consists in choosing s = p N/d rather than √ N , see §3.4.
3. Main results
In this section, we generalize the algorithms from §2, and show that they apply tothe general setting of q -holonomic sequences. 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 .11 efinition 1 ( q -holonomic sequence) . Let K be a field, and q ∈ K . A sequence ( u n ( q )) n ≥ in K N is called q -holonomic if there exist r ∈ N and polynomials c ( 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 ( q, q n ) u n ( q ) = 0 , for all n ≥ . (6) The integer r is called the order of the q -recurrence (6) . When r = 1 , we say that ( u n ( q )) n ≥ is q -hypergeometric. The most basic examples are the q -bracket and the q -factorial , [ n ] q := 1 + q + · · · + q n − and [ n ] q ! := n Y k =1 [ k ] q . (7)They are clearly q -holonomic, and even q -hypergeometric.The sequences ( u n ) n ≥ = ( q n ) n ≥ , ( v n ) n ≥ = ( q n ) n ≥ and ( w n ) n ≥ = ( q ( n )) n ≥ arealso q -hypergeometric, since they satisfy the recurrence relations u n +1 − qu n = 0 , v n +1 − q n +1 v n = 0 , w n +1 − q n w n = 0 . However, the sequence ( q n ) n ≥ is not q -holonomic [72, Ex. 2.2(b)]. More generally, thisalso holds for the sequence ( q n s ) n ≥ , for any s > , see [26, Th. 4.1] and also [66, Th. 1.1].Another basic example is the q -Pochhammer symbol ( x ; q ) n := n − Y 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, thesequence ( q ; q ) n := Q nk =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 = P nk =0 p k q ak + 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 ( x, y ) , . . . , c r ( x, y ) are constant in the variable y . Forthese kinds of sequences there exist quasi-optimal algorithms [41, 61, 107], therefore weassume from now on that the maximal degree d of c ( 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 (framedoriented) knot (in 3-space) is a powerful knot invariant, related to the Alexander poly-nomial [19]; it is a q -holonomic sequence of Laurent polynomials [71]. Its recurrenceequations are themselves of interest, as they are closely related to the A-polynomial of aknot, via the AJ conjecture [54, 65, 67], verified in some cases using massive computeralgebra calculations [69].It is well known that the class of q -holonomic sequences is closed under several op-erations, such as addition, multiplication, Hadamard product and monomial substitu-tion [72, 91, 96]. All these closure properties are effective, i.e. , they can be executedalgorithmically on the level of q -recurrences. Several computer algebra packages areavailable for the manipulation of q -holonomic sequences, e.g. , the Mathematica pack-ages qGeneratingFunctions [91] and HolonomicFunctions [98], and the Maple packages qsum [28], qFPS [132], qseries and
QDifferenceEquations .12 simple but useful fact is that the order- r scalar q -recurrence (6) can be translatedinto a first-order recurrence on r × vectors: u n + r ... u n +1 = − c r − c r · · · − c c r − c c r · · · ... . . . ... ... · · · × u n + r − ... u n . (9)In particular, the N -th term of the q -holonomic sequence ( u n ) is simply expressiblein terms of the matrix q -factorial M ( q N − ) · · · M ( q ) M (1) , (10)where M ( q n ) denotes the companion matrix from equation (9). This observation iscrucial, since it exposes the connection to the algorithms presented in the previous section. 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 coveredby [38, §6], where the cost O ( M ( √ N )) is reached under some invertibility assumptions. Theorem 1.
Let K be a field, let q ∈ K \ { } and N ∈ N . The q -factorial [ N ] q ! can becomputed using O ( M ( √ N )) operations in K . The same is true for the q -Pochhammersymbol ( α ; q ) N for any α ∈ K .Proof. If α = 0 , then ( α ; q ) N = 1 . If q = 0 , then [ N ] q ! = 1 and ( α ; q ) N = 1 − α . Wecan assume that q ∈ K \ { , } and α ∈ K \ { } . We have [ N ] q ! = r N · F ( q − ) and ( α ; q ) N = α N · F ( α − ) , where r := q/ (1 − q ) and F ( x ) := Q N − i =0 ( x − q i ) . Algorithm 2 can be used to compute F ( q − ) and F ( α − ) in O ( M ( √ N )) operations in K . The costof 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 computein O ( M ( √ N )) operations in K : • the q -binomial coefficient (cid:0) Nn (cid:1) q ; • the coefficient of x n in the polynomial Q Nk =1 (1 + q k − x ) ; • the sum (cid:0) N − n (cid:1) q (cid:0) n (cid:1) q + q (cid:0) N − n (cid:1) q (cid:0) n (cid:1) q + · · · + q n (cid:0) N − nn (cid:1) q (cid:0) nn (cid:1) q .Proof. The first assertion is a direct consequence of Theorem 1. The second assertion isa consequence of the first one, and of (1). The third assertion is a consequence of thefirst one, and of (2). 13 .3. N -th term of a q -holonomic sequence We now offer the promised q -analogue of Chudnovskys’ result on the computationof the N -th term of an arbitrary holonomic sequence in O ( M ( √ N ) log N ) arithmeticoperations. Note that Chudnovskys’ case q = 1 is also covered by [38, §6], where theimproved cost O ( M ( √ N )) is reached under additional invertibility assumptions. Theorem 3.
Let K be a field, q ∈ K \ { } and N ∈ N . Let ( u n ( q )) n ≥ be a q -holonomicsequence satisfying recurrence (6) , and assume that c r ( q, q k ) is nonzero for k = 0 , . . . , N − . 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 − ) · · · M ( q ) M (1) can be computed in O ( M ( √ N )) , where M ( q n ) denotes the companion matrix fromequation (9). 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 − ) · · · M ( q ) M (1) can be simplified using U N ( q ) = M ( q k ) · · · M (1) · U n ( q ) r , where r = ⌊ ( N − /n ⌋ and k = N − − 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 ofjust O ( M ( √ n ) + log( N )) arithmetic operations. It follows that if, for instance, the basefield K is the prime field F p , then the prime number p should be larger than N in orderto 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 ) := X n ≥ x n [ n ] q ! , and let E ( N ) q ( x ) := e q ( x ) mod x N be its polynomial truncation of degree N − . 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 ≥ is q -hypergeometric, andsatisfies the recurrence [ n +1] q u n +1 ( q ) − αu n = 0 , therefore v N ( q ) := P N − i =0 u i ( q ) satisfiesthe second-order recurrence [ n +1] q ( v n +2 ( q ) − v n +1 ( q )) − α ( v n +1 ( q ) − v n ( q )) = 0 . ApplyingTheorem 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 q -difference equation. For instance, one can evaluate fast all truncations of Heine’s q -hypergeometric series φ ([ a, b ] , [ c ]; q ; x ) := X n ≥ ( a ; q ) n ( b ; q ) n ( c ; q ) n · x n ( q ) n . .4. 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 complexityanalysis 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 matrixof degree d ≥ . We wish to examine the amount of field operations in K needed for thecomputation of M ( q N − ) · · · M (1) in terms of N, d and n . Recall that MM ( n, d ) controlsthe arithmetic complexity of the product in M n ( K [ x ] d ) and it holds that MM ( n, d ) = O ( n θ d + n M ( d )) = ˜ O ( 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 (1), where we compute P s ( x ) = M ( q s − x ) · · · M ( x ) using the divide-and-conquer scheme P t ( x ) = P t ( q t x ) · P t ( x ) and P t +1 ( x ) = M ( q t x ) · P t ( q t x ) · P t ( x ) . Note that P t ( x ) is an n × n polynomial matrix ofdegree at most td and therefore the cost of this step is O ( MM ( n, sd )) = ˜ O ( n θ d √ N ) .Step (2) is done component-wisely at each entry of P ( x ) . By the explained fast chirptransform algorithms, it essentially boils down to n multiplications of two polynomials,one of degree sd and the other of degree sd . The cost of the second step is therefore O ( n M ( sd )) = ˜ O ( n d √ N ) . The last step is the multiplication of N/s = s matrices withentries in K and has therefore an arithmetic complexity of O ( n θ s ) = ˜ O ( 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 = p N/d , which we again can assumeto be integer . Then, by the same arguments as above, the costs of the three steps are O ( MM ( n, √ N d )) = ˜ O ( n θ √ N d ) , O ( n M ( √ N d )) = ˜ O ( n √ N d ) and O ( n θ √ N d ) respec-tively.Now, we address specifically the computation of the N -th term in a q -holonomicsequence. If ( u n ( q )) n ≥ is given by a q -recurrence c r ( q, q n ) u n + r ( q ) + · · · + c ( 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 c r ( q, q N − ) · · · c r ( q, q ) c r ( q, · (cid:2) · · · (cid:3) × ˜ M ( q N − ) · · · ˜ M ( q ) ˜ M (1) × ... , where now ˜ M ( x ) := c r ( q, x ) · M ( x ) = − c r − ( q, x ) · · · − c ( q, x ) − c ( q, x ) c r ( q, x ) · · · . . . ... · · · c r ( q, x ) 0 . Similarly as before, if p N/d is not an integer, then we can compute u N ( q ) first, where N = ⌊ p N/d ⌋ d , and then proceed “naively”. Note that N − N < √ Nd − d = O ( √ Nd ) . ˜ M ( q N − ) · · · ˜ M (1) and c r ( q, q N − ) · · · c r ( q, . If the degreesof c ( q, y ) , . . . , c r ( q, y ) are bounded by d , then the considerations above imply that thetwo q -factorials can be computed in O ( MM ( r, √ N d )+ r M ( √ N d )) and O ( M ( √ N d )) op-erations in K , respectively. We obtain the following theorem (compare with [37, Thm. 2]). Theorem 5.
Under the assumptions of Theorem 3, let d ≥ be the maximum of thedegrees of c ( q, y ) , . . . , c r ( q, y ) . Then, for any N > d , the term u N ( q ) can be computedin O ( r θ √ N d + r M ( √ N d )) operations in K . Theorem 3 can be adapted to the computation of several coefficients of a q -holonomicsequence. The proof is similar to that of Theorem 15 in [38], however simpler, becausewe deal with geometric progressions instead of arithmetic ones. Theorem 6.
Under the assumptions of Theorem 5, let N < N < · · · < N n = N bepositive integers, where n ≤ √ N . Then, the terms u N ( q ) , . . . , u N n ( q ) can be computedaltogether in O ( M ( √ N ) log N ) operations in K .Proof. As before, we assume that N is a perfect square; let s = √ N . Examining thepresented algorithms, we notice that on the way of computing the matrix q -factorial U N := M ( q N − ) · · · M ( q ) M (1) we obtain the evaluated polynomials P (1) , P ( Q ) , . . . , P ( Q s − ) , where Q := q s and P ( x ) := M ( q s − x ) · · · M ( qx ) M ( x ) . The q -factorial U N is then foundby step (3) by trivially multiplying P ( Q s − ) · · · P ( Q ) P (1) . Observe that while multiply-ing together from right to left we actually also automatically compute P ( Q j − ) · · · P ( Q ) P (1) = M ( q sj − ) · · · M ( q ) M (1) = U sj , for every j = 1 , . . . , s . It follows that employing Algorithm 3 and by simply taking thetop right element of each U sj , we find not only u N , but actually u s , u s , . . . , u s = u N .This already indicates that simultaneous computation of s terms is achievable in similarcomplexity after some “distillation”. In general, we are interested in the sequence of u i ( q ) at indices i = N , . . . , N n , hence we need to perform the following refinement step.Let d ∈ N be a positive integer with d ≤ n ≤ √ N and assume that for some k (0)1 , . . . , k (0) n with k (0) j ≤ N j < k (0) j + 2 d 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 polynomialmatrix P d ( x ) = M ( q d − x ) · · · M ( qx ) M ( x ) . Compute then the values q k (0)1 , . . . , q k (0) n andevaluate P d ( x ) simultaneously at them. For each j = 1 , . . . , n it holds that P d ( q k (0) j ) · U k (0) j = U k (0) j + d . We perform this multiplication for those indices j for which k (0) j + d ≤ N j < k (0) j + 2 d .For these j we then set k (1) j = k (0) j + d and let k (1) j = k (0) j for the other indices; moreover d := ⌈ d / ⌉ . We iterate this process at most ℓ := ⌈ log( d ) ⌉ many times until d ℓ = 1 .Then we can easily find U N , . . . , U N n , from which we finally deduce u N , . . . , u N n .16ach 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 , U s , . . . , U s in O ( M ( √ N )) base operations (Theorem 5). Hence, we may choose d = s and for each j = 1 , . . . , n let k (0) j be the largest element in { , s, . . . , ( s − s } suchthat k (0) j ≤ N j . Clearly, all conditions of the above refinement step are satisfied and weneed at most ⌈ log( s ) ⌉ = O (log N ) many such steps. The total complexity is henceforth 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 ε forsome ε > , then we can omit the log -factor in N : Corollary 7.
Under the assumptions of Theorem 5, let N < N < · · · < N n = N be pos-itive integers, where n < N − ε for some < ε < . Then, the terms u N ( q ) , . . . , u N n ( 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, regardingcomplexity, we use O ( M ( n )) = O ( M ( N − ε )) = O ( M ( √ N ) N − ε ) and obtain that thesame 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 severalcoefficients, we have the following trade-off: either we compute at most √ N terms in thearithmetic complexity O (( r θ d √ N + r M ( d √ N )) log N ) , or at most p N/d terms, but ina better cost of O (( r θ √ N d + r M ( √ N d )) log N ) . The proofs combine the considerationsabove with setting s = √ N and s = p N/d respectively. In both cases we can get rid ofthe log-factor in N like in Corollary 7 by computing a factor of N ε less terms. q is an integer: bit complexity Until now, we only considered the arithmetic complexity model, which is very well-suited to measure the algorithmic cost when working in algebraic structures whose basicinternal 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 arith-metic complexity model needs to be replaced by the bit-complexity model.Recall that the most basic operations on integer numbers can be performed in quasi-optimal time, that is, in a number of bit operations which is almost linear, up to loga-rithmic factors, in (the maximum of) their bit size. The most basic operation is integermultiplication, for which quasi-linear time algorithms are known since the early seven-ties, starting with the famous paper by Schönhage and Strassen [130] who showed thattwo n -bit integers can be multiplied in O ( n log n log log n ) bit operations. After severalsuccessive improvements, e.g., [62, 81], we know as of 2020 that two n -bit integers can bemultiplied in time O ( n log n ) [80]. 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 inthe usual holonomic setting. For example, the bitsize of u N = N ! is O ( N log N ) , how-ever both, the “naive” method of computing it using u n = nu n − , or Strassen’s baby-steps/giant-steps method, yield worse bit complexity. In the naive approach the problem17s that integers of unbalanced bitsize are multiplied together and hence not the full powerof fast integer multiplication techniques can be employed. A more clever and very simpleway is to just use the fact that N ! = (1 · · · ⌊ N/ ⌋ ) × (( ⌊ N/ ⌋ + 1) · · · N ) , and that the bitsizes of both factors have magnitude O ( N/ N )) = ˜ O ( N ) . Thus,fast integer multiplication can be used to multiply them in ˜ O ( N ) bit-complexity. Thisidea results in the binary splitting Algorithm 5 . This algorithm is very classical, and weonly recall it for completeness. See [23, §12] for a good survey on this technique, and itsapplications.
Algorithm 5 ( BinSplit ) Input : A = [ a , . . . , a N ] list of elements from some arbitrary ring R Output : a N · · · a function F ( A ) if N = 1 then return A [1] return F ( A [ ⌊ N/ ⌋ + 1 , . . . , N ]) · F ( A [1 , . . . , ⌊ N/ ⌋ ]) Assume that each a i in the input of BinSplit has at most k bits and let C ( n ) be thecomplexity of BinSplit if A = [ a , . . . , a n ] is n -dimensional. It follows that C ( N ) ≤ C ( ⌈ N/ ⌉ ) + M Z ( N/ · k ) , where M Z ( n ) = ˜ O ( n ) is the cost of multiplication of integers with n bits. We obtain C ( N ) = ˜ O ( N k ) . Hence, using this method, the computation of u N = N ! has ˜ O ( N ) bit-complexity, which is quasi-optimal. Moreover, the same idea applies to any holo-nomic sequence, by deducing the first order matrix recurrence and computing the matrixproduct using BinSplit .Now we shall see that the q -holonomic case is similar. First, let q be a positive integerof B bits and consider the computation of the q -factorial u N ( q ) = (1 + q )(1 + q + q ) · · · (1 + q + · · · + q N − ) , as an illustrative example. For each factor we have the trivial inequalities q n < q + · · · + q n < q n +1 , meaning that q N ( N − / < u N ( q ) < q N ( N +1) / , so the bitsize of u N ( q ) is of magni-tude N B . The “naive” algorithm of deducing u N ( q ) by first computing the integers q i , then the corresponding sums and products, has ˜ O ( N B ) binary complexity. Thismethod is not (quasi-)optimal with respect to the output size. It is also easy to see thatthe presented baby-steps/giant-steps based algorithms yield bad bit-complexity as well,despite their good arithmetic cost.Similarly, if u N ( q ) = N − X n =0 q n , u N ( q ) is bounded in absolute value from above by N q ( N − and by q ( N − from below, so its bitsize is again of magnitude N B . The “naive” algorithmconsisting of computing the terms q i one after the other before summing, has againnon-optimal bit-complexity ˜ O ( N B ) .Can one do better? The answer is “yes” and one can even achieve a complexitywhich is quasi-linear in the bitsize of the output. Similarly to the holonomic setting, itis sufficient to use the q -holonomic character of u N ( q ) , and to reduce its computationto that of a q -factorial matrix as in §2.2, which can then be handled with BinSplit . Tobe more precise, given an integer or rational number q of bitsize B and any q -holonomicsequence ( u n ( q )) n ≥ defined by polynomials c , . . . , c r ∈ Z [ x, y ] with c r ( q, q n ) = 0 forany n ∈ N , we define M ( x ) ∈ Q ( x ) by − c r − ( q,x ) c r ( q,x ) · · · − c ( q,x ) c r ( q,x ) − c ( q,x ) c r ( q,x ) · · · ... . . . ... ... · · · . Then, as observed before, u N can be read off from M ( q N − ) · · · M ( q ) M (1) , which we aim to compute efficiently. Again, instead of using baby-steps/giant-steps, it isa better idea to use binary splitting by applying Algorithm 5 to A = [ M (1) , . . . , M ( q N − )] .Note that obviously, any element in A is a matrix with rational entries of bitsize boundedby O ( N B ) . Therefore, the complexity of BinSplit does not exceed ˜ O ( N B ) by the sameargument as before, now using fast multiplication of rational numbers. These consider-ations prove Theorem 8.
Under the assumptions of Theorem 3, with K = Q , the term u N ( q ) can becomputed in ˜ O ( N B ) bit operations, where B is the bitsize of q . As a corollary, (truncated) solutions of q -difference equations can be evaluated usingthe same (quasi-linear) bit-complexity. This result should be viewed as the q -analogue ofthe 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 (withoutproof) by Schroeppel and Salamin in Item 178 of [21].
4. Applications q -holonomic sequences As already mentioned, many q -holonomic sequences arise in combinatorics, for ex-ample in connection with the enumeration of lattice polygons, where q -analogues of theCatalan numbers n +1 (cid:0) nn (cid:1) occur naturally [63, 76], or in the enumeration of special fami-lies of matrices with coefficients in the finite field F q [93, 94, 142], where sequences relatedto the Gaussian coefficients (cid:0) nk (cid:1) q also show up.19 huge subfield of combinatorics is the theory of partitions [14], where q -holonomicsequences occur as early as in the famous Rogers-Ramanujan identities [123, 124], seealso [14, Ch. 7], e.g. , X n ≥ q n (1 − q ) · · · (1 − q n ) = Y n ≥ − q n +1 )(1 − q n +4 ) which translates the fact that the number of partitions of n into parts that differ by atleast is equal to the number of partitions of n into parts congruent to or modulo . Andrews [12, 13], see also [14, Chapter 8], laid the foundations of a theory able tocapture 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 comingfrom partitions can be evaluated fast. For instance, the family of truncated polynomials F n ( x ) := ∞ Y k =1 (1 − x k ) mod x n , can be evaluated fast due to our results and to the identity [113, §6] F N ( q ) = X ( n +12 ) 5. Experiments Algorithms 1 and were implemented in Magma and Algorithm 3 in Maple . All im-plementations deliver some encouraging timings. Of course, since these algorithms aredesigned to be fast in the arithmetic model , it is natural to make experiments over afinite field K , or over truncations of real/complex numbers, as was done in [110] for theproblem in §2.2.Recall that both Algorithms 1 and compute Q N − i =0 ( α − q i ) ∈ K , given α, q ina field K , and N ∈ N , whereas Algorithm 3 finds M ( q N − ) · · · M ( q ) M (1) ∈ K n × n for a23egree N Naive algorithm Algorithm 1 Algorithm 2 . 04 0 . 03 0 . . 18 0 . 03 0 . . 72 0 . 06 0 . . 97 0 . 14 0 . . 79 0 . 32 0 . . 16 0 . 73 0 . . 56 1 . 68 0 . . 65 3 . 84 0 . . 25 8 . 65 0 . . 65 1 . . 42 2 . . 27 6 . . 58 14 . . 03 29 . . 51 61 . . 28 137 . . . . . Table 1 Comparative timings (in seconds) for the computation of Q N − i =0 ( α − q i ) ∈ F p , with p = 2 + 3 and ( α, q ) randomly chosen in F p × F p . All algorithms were executed on the samemachine, running Magma v. 2.24 . For each target degree N , each execution was limited to 1 hour. Naive algorithm could reach degree N = 2 , Algorithm 1 degree N = 2 , and Algorithm 2 degree N = 2 = 8 014 398 509 481 984 . By extrapolation, the Naive algorithm would have needed ≈ × . sec. ≈ years on the same instance, and Algorithm 2 approximately hours. 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 + 3 elements.Timings for Algorithms 1 and are presented in Table 1. We compare the straight-forward 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 onmultipoint 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 increasedfrom k to k +2 , timings are also multiplied (roughly) by 4 in column Naive , and(roughly) by 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 , ≈ of the time isspent in step (3) (resultant computation), the other steps taking ≈ each; for Algorithm 2 , step (1) takes ≈ , step (2) takes ≈ , and step (3) is negligible.24 − − − N T i m e i n s ec o nd s d = 1; n = 2 ; naive d = 1 ; n = 2 d = 1 ; n = 4 d = 1 ; n = 8 Figure 1: Timings of Algorithm 3 implemented in Maple 2020.2 . We compare d = 1 and n = 2 , , forvarious values of N . In order to visualize the performance of Algorithm 3 , we took random polynomialmatrices in M n ( F p [ x ]) of degree d . Figure 1 compares the time needed to compute theMatrix q -factorial for d = 1 and n = 2 , , as N growths. The black lines represent thebest linear fits to the data points and are given by the equations y = 0 . x − . , y =0 . x − . and y = 0 . x − . respectively. Note that we are plotting on a log-logscale, therefore the established complexity of ˜ O ( N / ) is indicated by the coefficient of x in these linear fits which is always only slightly greater than / . In the same figure, wealso show the timings for d = 1 and n = 2 of the naive algorithm given by successivelycomputing and multiplying M ( q i ) together. The best linear fit almost perfectly describesthis data and has a slope of . ≈ , in line with the linear complexity in N .In Figure 2 we show similar timings, however now for n = 2 and d = 2 , , . The linearfits to the data are now given by y = 0 . x − . , y = 0 . x − . and y = 0 . x − . .Again, they describe the observations very well and the coefficients of the regressions arein line with the proven complexity. We observe that, as expected, the lines have slopesof roughly / and are closer together than in the previous figure. The naive methodhas a slope of . ≈ . 6. 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 theholonomic (“ q = 1 ”) case. We have demonstrated through several examples that thisbasic algorithmic improvement has many other algorithmic implications, notably on thefaster evaluation of many families of polynomials and on the acceleration of algorithmsfor q -difference equations.Here are some questions that should be investigated in the future.25 − − − N T i m e i n s ec o nd s d = n = 2 ; naive d = 2 ; n = 2 d = 4 ; n = 2 d = 8 ; n = 2 Figure 2: Timings of Algorithm 3 implemented in Maple 2020.2 . We compare n = 2 and d = 2 , , forvarious values of N . 1. (Counting points on q -curves) Counting efficiently points on (hyper-)elliptic curvesleads to questions like: for a, b ∈ Z , compute the coefficient of x p − in G p ( x ) :=( x + ax + b ) p − modulo p , for one [38] or several [79] primes p . A natural extensionis to ask the same with G p ( x ) replaced by Q p − k =1 ( q k x + aq k x + b ) . This mighthave applications related to §4.4, or to counting points on q -deformations [128].2. (Computing q -deformed real numbers) Recently, Morier-Genoud and Ovsienko [108]introduced q -analogues of real numbers, see also [103, 109]. How fast can one com-pute (truncations or evaluations of) quantized versions of numbers like e or π ?3. (Evaluating more polynomials) Is it possible to evaluate fast polynomials of theform P Nℓ =0 x ℓ s , for s ≥ , and many others that escape the q -holonomic class? E.g. , [24] presents a beautiful generalization of Algorithm 1 to the fast evaluationof isogenies between elliptic curves, by using elliptic resultants , with applicationsto isogeny-based cryptography.4. (“Precise” complexity of q -holonomic sequences) Can one prove non-trivial lowerbounds, ideally matching the upper bounds, on examples treated in this article? Acknowledgements. We would like to thank Luca De Feo for his initial question, whomotivated this work, and for the very interesting subsequent discussions. Our warmthanks go to Lucia Di Vizio for valuable comments that led to §4.4. Moreover, we thankher 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 thesecond author by the Austrian Science Fund (FWF) P-31338.26 eferences [1] S. Abramov, M. Bronstein, M. Petkovšek, On polynomial solutions of linear operator equations,in: ISSAC’95, ACM Press, 1995, pp. 290–296.[2] S. A. Abramov, Rational solutions of linear difference and q -difference equations with polynomialcoefficients, Programmirovanie (6) (1995) 3–11.[3] S. A. Abramov, A direct algorithm to compute rational solutions of first order linear q -differencesystems, Discrete Math. 246 (1-3) (2002) 3–12.URL https://doi.org/10.1016/S0012-365X(01)00248-5 [4] S. A. Abramov, P. Paule, M. Petkovšek, q -hypergeometric solutions of q -difference equations,Discrete Math. 180 (1-3) (1998) 3–22.URL https://doi.org/10.1016/S0012-365X(97)00106-4 [5] S. A. Abramov, E. V. Zima, D’Alembertian solutions of inhomogeneous linear equations (differ-ential, difference, and some other), in: ISSAC’96, ACM Press, 1996, p. 232–240.URL https://doi.org/10.1145/236869.237080 [6] B. Adamczewski, J. P. Bell, E. Delaygue, F. Jouhet, Congruences modulo cyclotomic polynomialsand algebraic independence for q -series, Sém. Lothar. Combin. 78B (2017) Art. 54, 12.[7] C. R. Adams, On the linear ordinary q -difference equation, Ann. of Math. (2) 30 (1-4) (1928/29)195–205.URL https://doi.org/10.2307/1968274 [8] C. R. Adams, Linear q -difference equations, Bull. AMS 37 (6) (1931) 361–400.URL https://doi.org/10.1090/S0002-9904-1931-05162-4 [9] W. A. Al-Salam, L. Carlitz, Some orthogonal q -polynomials, Math. Nachr. 30 (1965) 47–61.URL https://doi.org/10.1002/mana.19650300105 [10] M. Aldaz, G. Matera, J. L. Montaña, L. M. Pardo, A new method to obtain lower bounds forpolynomial evaluation, Theoret. Comput. Sci. 259 (1-2) (2001) 577–596.URL https://doi.org/10.1016/S0304-3975(00)00149-3 [11] J. Alman, V. Vassilevska Williams, A refined laser method and faster matrix multiplication, arXivpreprint, 29 pages (2020).URL https://arxiv.org/abs/2010.05846 [12] G. E. Andrews, Partition identities, Advances in Math. 9 (1972) 10–51.URL https://doi.org/10.1016/0001-8708(72)90028-X [13] G. E. Andrews, A general theory of identities of the Rogers-Ramanujan type, Bull. AMS 80 (1974)1033–1052.URL https://doi.org/10.1090/S0002-9904-1974-13616-5 [14] G. E. Andrews, The theory of partitions, Addison-Wesley Publishing Co., Reading„ 1976, ency-clopedia of Mathematics and its Applications, Vol. 2.[15] G. E. Andrews, The fifth and seventh order mock theta functions, Trans. Amer. Math. Soc. 293 (1)(1986) 113–134.URL https://doi.org/10.2307/2000275 [16] G. E. Andrews, q -Catalan identities, in: The legacy of Alladi Ramakrishnan in the mathematicalsciences, Springer, 2010, pp. 183–190.URL https://doi.org/10.1007/978-1-4419-6263-8_10 [17] C. E. Arreche, Y. Zhang, Computing differential Galois groups of second-order linear q -differenceequations, preprint, 2020.[18] R. Askey, Continuous q -Hermite polynomials when q > , in: q -series and partitions, vol. 18 ofIMA Vol. Math. Appl., Springer, 1989, pp. 151–158.URL https://doi.org/10.1007/978-1-4684-0637-5_12 [19] D. Bar-Natan, S. Garoufalidis, On the Melvin-Morton-Rozansky conjecture, Invent. Math. 125 (1)(1996) 103–133.URL https://doi.org/10.1007/s002220050070 [20] M. Barkatou, T. Cluzeau, A. El Hajj, Simple forms and rational solutions of pseudo-linear systems,in: ISSAC’19, ACM Press, 2019, pp. 26–33.URL https://doi.org/10.1145/3326229.3326246 https://doi.org/10.1017/s0025557200044491 23] D. J. Bernstein, Fast multiplication and its applications, in: Algorithmic number theory: lattices,number fields, curves and cryptography, vol. 44 of Math. Sci. Res. Inst. Publ., Cambridge Univ.Press, 2008, pp. 325–384.[24] D. J. Bernstein, L. D. Feo, A. Leroux, B. Smith, Faster computation of isogenies of large primedegree, Cryptology ePrint Archive, Report 2020/341, https://eprint.iacr.org/2020/341 , pre-sented to ANTS-XIV (2020).[25] J.-P. Bézivin, Les suites q -récurrentes linéaires, Compositio Math. 80 (3) (1991) 285–307.URL [26] J.-P. Bézivin, Sur les équations fonctionnelles aux q -différences, Aequationes Math. 43 (2-3) (1992)159–176.URL https://doi.org/10.1007/BF01835698 [27] L. I. Bluestein, A linear filtering approach to the computation of the discrete Fourier transform,IEEE Trans. Electroacoustics AU-18 (1970) 451–455.[28] H. Böing, W. Koepf, Algorithms for q -hypergeometric summation in computer algebra, J. SymbolicComput. 28 (6) (1999) 777–799.URL https://doi.org/10.1006/jsco.1998.0339 [29] A. Borodin, S. Cook, On the number of additions to compute specific polynomials, SIAM J.Comput. 5 (1) (1976) 146–157.URL https://doi.org/10.1137/0205013 [30] P. B. Borwein, Padé approximants for the q -elementary functions, Constr. Approx. 4 (4) (1988)391–402.URL https://doi.org/10.1007/BF02075469 [31] A. Bostan, Computing the N -th Term of a q -Holonomic Sequence, in: ISSAC’20, ACM Press,2020, pp. 46–53.[32] A. Bostan, X. Caruso, E. Schost, A fast algorithm for computing the characteristic polynomial ofthe p -curvature, in: ISSAC’14, ACM Press, 2014, pp. 59–66.URL https://doi.org/10.1145/2608628.2608650 [33] A. Bostan, X. Caruso, E. Schost, A fast algorithm for computing the p -curvature, in: ISSAC’15,ACM, New York, 2015, pp. 69–76.[34] A. Bostan, X. Caruso, E. Schost, Computation of the similarity class of the p -curvature, in:ISSAC’16, ACM Press, 2016, pp. 111–118.[35] A. Bostan, F. Chyzak, T. Cluzeau, B. Salvy, Low complexity algorithms for linear recurrences, in:ISSAC’06, ACM Press, 2006, pp. 31–38.URL https://doi.org/10.1145/1145768.1145781 [36] A. Bostan, F. Chyzak, M. Giusti, R. Lebreton, G. Lecerf, B. Salvy, E. Schost, Algorithmes Efficacesen Calcul Formel, Palaiseau, 2017, 686 pages, in French. Printed by CreateSpace. Also availablein electronic form.URL https://hal.archives-ouvertes.fr/AECF/ [37] A. Bostan, T. Cluzeau, B. Salvy, Fast algorithms for polynomial solutions of linear differentialequations, in: ISSAC’05, ACM Press, 2005, pp. 45–52.URL https://doi.org/10.1145/1073884.1073893 [38] A. Bostan, P. Gaudry, E. Schost, Linear recurrences with polynomial coefficients and applicationto integer factorization and Cartier-Manin operator, SIAM J. Comput. 36 (6) (2007) 1777–1806.URL https://doi.org/10.1137/S0097539704443793 [39] A. Bostan, M. Kauers, Automatic classification of restricted lattice walks, in: FPSAC 2009, Dis-crete Math. Theor. Comput. Sci. Proc., AK, Assoc. Discrete Math. Theor. Comput. Sci., Nancy,2009, pp. 201–215.[40] A. Bostan, G. Lecerf, É. Schost, Tellegen’s principle into practice, in: Proceedings of ISSAC’03,ACM Press, 2003, pp. 37–44.[41] A. Bostan, R. Mori, A Simple and Fast Algorithm for Computing the N -th Term of a LinearlyRecurrent Sequence, in: SOSA’21 (Symposium on Simplicity in Algorithms), SIAM, 2021.URL https://specfun.inria.fr/bostan/BoMo21.pdf [42] A. Bostan, E. Schost, Polynomial evaluation and interpolation on special sets of points, J. Com-plexity 21 (4) (2005) 420–446.URL https://doi.org/10.1016/j.jco.2004.09.009 [43] A. Bostan, E. Schost, Fast algorithms for differential equations in positive characteristic, in: IS-SAC’09, ACM Press, 2009, pp. 47–54.URL https://doi.org/10.1145/1576702.1576712 [44] M. Bousquet-Mélou, Convex polyominoes and algebraic languages, J. Phys. A 25 (7) (1992) 1935– http://stacks.iop.org/0305-4470/25/1935 [45] P. Bürgisser, M. Clausen, M. A. Shokrollahi, Algebraic complexity theory, vol. 315 of Grundlehrender Mathematischen Wissenschaften, Springer, 1997.URL https://doi.org/10.1007/978-3-662-03338-8 [46] D. G. Cantor, E. Kaltofen, On fast multiplication of polynomials over arbitrary algebras, ActaInform. 28 (7) (1991) 693–701.[47] R. D. Carmichael, The General Theory of Linear q -Difference Equations, Amer. J. Math. 34 (2)(1912) 147–168.URL https://doi.org/10.2307/2369887 [48] P. Cartier, Démonstration “automatique” d’identités et fonctions hypergéométriques (d’après D.Zeilberger), Astérisque (206) (1992) Exp. No. 746, 3, 41–91, séminaire Bourbaki, Vol. 1991/92.[49] D. V. Chudnovsky, G. V. Chudnovsky, Approximations and complex multiplication according toRamanujan, in: Ramanujan revisited (Urbana-Champaign, Ill., 1987), Academic Press, Boston,MA, 1988, pp. 375–472.[50] F. Chyzak, Gröbner bases, symbolic summation and symbolic integration, in: Gröbner bases andapplications (Linz, 1998), vol. 251 of London Math. Soc. Lecture Note Ser., Cambridge Univ.Press, 1998, pp. 32–60.[51] F. Chyzak, An extension of Zeilberger’s fast algorithm to general holonomic functions, DiscreteMath. 217 (1-3) (2000) 115–134.URL https://doi.org/10.1016/S0012-365X(99)00259-9 [52] F. Chyzak, P. Dumas, H. Le, J. Martin, M. M., B. Salvy, Taming apparent singularities via Oreclosure, manuscript, 2016.[53] E. Costa, R. Gerbicz, D. Harvey, A search for Wilson primes, Math. Comp. 83 (290) (2014) 3071–3091.URL https://doi.org/10.1090/S0025-5718-2014-02800-7 [54] R. Detcherry, S. Garoufalidis, A diagrammatic approach to the AJ conjecture, Math. Ann. 378 (1-2) (2020) 447–484.URL https://doi.org/10.1007/s00208-020-02028-y [55] L. Di Vizio, Arithmetic theory of q -difference equations: the q -analogue of Grothendieck-Katz’sconjecture on p -curvatures, Invent. Math. 150 (3) (2002) 517–578.URL https://doi.org/10.1007/s00222-002-0241-z [56] L. Di Vizio, C. Hardouin, Intrinsic approach to Galois theory of q -difference equations, with thepreface to Part IV by Anne Granier, Mem. Amer. Math. Soc. 77 pages. To appear.URL http://arxiv.org/abs/1002.4839 [57] L. Di Vizio, J.-P. Ramis, J. Sauloy, C. Zhang, Équations aux q -différences, Gaz. Math. (96) (2003)20–49.[58] T. Ekedahl, G. van der Geer, Cycle classes on the moduli of K3 surfaces in positive characteristic,Selecta Math. (N.S.) 21 (1) (2015) 245–291.URL https://doi.org/10.1007/s00029-014-0156-8 [59] S. B. Ekhad, D. Zeilberger, The number of solutions of X = 0 in triangular matrices over GF( q ) ,Electron. J. Combin. 3 (1) (1996) Research Paper 2, approx. 2.URL [60] T. Ernst, A comprehensive treatment of q -calculus, Birkhäuser/Springer, 2012.URL https://doi.org/10.1007/978-3-0348-0431-8 [61] C. M. Fiduccia, An efficient formula for linear recurrences, SIAM J. Comput. 14 (1) (1985) 106–112.URL https://doi.org/10.1137/0214007 [62] M. Fürer, Faster integer multiplication, SIAM J. Comput. 39 (3) (2009) 979–1005.URL https://doi.org/10.1137/070711761 [63] J. Fürlinger, J. Hofbauer, q -Catalan numbers, J. Combin. Theory Ser. A 40 (2) (1985) 248–264.URL https://doi.org/10.1016/0097-3165(85)90089-5 [64] F. L. Gall, Powers of tensors and fast matrix multiplication, in: ISSAC’14, 2014, pp. 296–303.[65] S. Garoufalidis, On the characteristic and deformation varieties of a knot, in: Proceedings of theCasson Fest, vol. 7 of Geom. Topol. Monogr., Geom. Topol. Publ., Coventry, 2004, pp. 291–309.URL https://doi.org/10.2140/gtm.2004.7.291 [66] S. Garoufalidis, The degree of a q -holonomic sequence is a quadratic quasi-polynomial, Electron.J. Combin. 18 (2) (2011) Paper 4, 23.URL https://doi.org/10.37236/2000 67] S. Garoufalidis, Quantum knot invariants, Res. Math. Sci. 5 (1) (2018) Paper No. 11, 17.URL https://doi.org/10.1007/s40687-018-0127-3 [68] S. Garoufalidis, C. Koutschan, Twisting q -holonomic sequences by complex roots of unity, in:ISSAC’12, ACM Press, 2012, pp. 179–186.URL https://doi.org/10.1145/2442829.2442857 [69] S. Garoufalidis, C. Koutschan, Irreducibility of q -difference operators and the knot , Algebr.Geom. Topol. 13 (6) (2013) 3261–3286.URL https://doi.org/10.2140/agt.2013.13.3261 [70] S. Garoufalidis, A. D. Lauda, T. T. Q. Lê, The colored HOMFLYPT function is q -holonomic,Duke Math. J. 167 (3) (2018) 397–447.URL https://doi.org/10.1215/00127094-2017-0030 [71] S. Garoufalidis, T. T. Q. Lê, The colored Jones function is q -holonomic, Geom. Topol. 9 (2005)1253–1293.URL https://doi.org/10.2140/gt.2005.9.1253 [72] S. Garoufalidis, T. T. Q. Lê, A survey of q -holonomic functions, Enseign. Math. 62 (3-4) (2016)501–525.URL https://doi.org/10.4171/LEM/62-3/4-7 [73] F. G. Garvan, New fifth and seventh order mock theta function identities, Ann. Comb. 23 (3-4)(2019) 765–783.URL https://doi.org/10.1007/s00026-019-00438-7 [74] J. von zur Gathen, J. Gerhard, Modern computer algebra, 3rd ed., Cambridge Univ. Press, 2013.URL https://doi.org/10.1017/CBO9781139856065 [75] C. F. Gauss, Summatio quarundam serierum singularium, Opera, Vol. 2, Göttingen: Gess. d. Wiss.(1863).[76] I. Gessel, A noncommutative generalization and q -analog of the Lagrange inversion formula, Trans.Amer. Math. Soc. 257 (2) (1980) 455–482.URL https://doi.org/10.2307/1998307 [77] W. Hahn, Über die höheren Heineschen Reihen und eine einheitliche Theorie der sogenanntenspeziellen Funktionen, Math. Nachr. 3 (1950) 257–294.URL https://doi.org/10.1002/mana.19490030502 [78] G. Hanrot, M. Quercia, P. Zimmermann, The middle product algorithm. I, Appl. Algebra Engrg.Comm. Comput. 14 (6) (2004) 415–438.URL https://doi.org/10.1007/s00200-003-0144-2 [79] D. Harvey, Counting points on hyperelliptic curves in average polynomial time, Ann. of Math. (2)179 (2) (2014) 783–803.URL https://doi.org/10.4007/annals.2014.179.2.7 [80] D. Harvey, J. van der Hoeven, Integer multiplication in time O ( n log n ) , Ann. of Math. To appear.[81] D. Harvey, J. van der Hoeven, G. Lecerf, Even faster integer multiplication, J. Complexity 36(2016) 1–30.URL https://doi.org/10.1016/j.jco.2016.03.001 [82] E. Heine, Untersuchungen über die Reihe (1 − q α )(1 − q β )(1 − q )(1 − q γ ) x + (1 − q α )(1 − q α +1 )(1 − q β )(1 − q β +1 )(1 − q )(1 − q )(1 − q γ )(1 − q γ +1 ) x + .. ,J. reine angew. Math. 34 (1847) 285–328.[83] J. Heintz, M. Sieveking, Lower bounds for polynomials with algebraic coefficients, Theoret. Com-put. Sci. 11 (3) (1980) 321–330.URL https://doi.org/10.1016/0304-3975(80)90019-5 [84] P. A. Hendriks, An algorithm for computing a standard form for second-order linear q -differenceequations, J. Pure Appl. Algebra 117/118 (1997) 331–352.URL https://doi.org/10.1016/S0022-4049(97)00017-0 [85] J. Hua, Counting representations of quivers over finite fields, J. Algebra 226 (2) (2000) 1011–1033.URL https://doi.org/10.1006/jabr.1999.8220 [86] M. E. H. Ismail, Lectures on q -orthogonal polynomials, in: Special functions 2000: current per-spective and future directions (Tempe, AZ), vol. 30 of NATO Sci. Ser. II Math. Phys. Chem.,Kluwer Acad. Publ., 2001, pp. 179–219.URL https://doi.org/10.1007/978-94-010-0818-1_8 [87] F. H. Jackson, On q -functions and a certain difference operator, Trans. Roy. Soc. Edin. 46 (11)(1909) 253–281.URL https://doi.org/10.1017/S0080456800002751 [88] F. H. Jackson, On q -definite integrals, Quar. J. Pure Appl. Math. 41 (1910) 193–203. 89] F. H. Jackson, q -Difference Equations, Amer. J. Math. 32 (4) (1910) 305–314.URL https://doi.org/10.2307/2370183 [90] N. M. Katz, Algebraic solutions of differential equations ( p -curvature and the Hodge filtration),Invent. Math. 18 (1972) 1–118.URL https://doi.org/10.1007/BF01389714 [91] M. Kauers, C. Koutschan, A Mathematica package for q -holonomic sequences and power series,Ramanujan J. 19 (2) (2009) 137–150.URL https://doi.org/10.1007/s11139-008-9132-2 [92] D. E. Khmel ′ nov, Improved algorithms for solving difference and q -difference equations, Program-mirovanie (2) (2000) 70–78.URL https://doi.org/10.1007/BF02759198 [93] A. A. Kirillov, On the number of solutions of the equation X = 0 in triangular matrices over afinite field, Funktsional. Anal. i Prilozhen. 29 (1) (1995) 82–87.URL https://doi.org/10.1007/BF01077044 [94] A. A. Kirillov, A. Melnikov, On a remarkable sequence of polynomials, in: Algèbre non commuta-tive, groupes quantiques et invariants (Reims, 1995), vol. 2 of Sémin. Congr., Soc. Math. France,Paris, 1997, pp. 35–42.[95] R. Koekoek, P. A. Lesky, R. F. Swarttouw, Hypergeometric orthogonal polynomials and their q -analogues, Monographs in Mathematics, Springer, 2010.URL https://doi.org/10.1007/978-3-642-05014-5 [96] W. Koepf, P. M. Rajković, S. D. Marinković, Properties of q -holonomic functions, J. DifferenceEqu. Appl. 13 (7) (2007) 621–638.URL https://doi.org/10.1080/10236190701264925 [97] T. H. Koornwinder, On Zeilberger’s algorithm and its q -analogue, J. Comput. Appl. Math. 48 (1-2)(1993) 91–111.URL https://doi.org/10.1016/0377-0427(93)90317-5 [98] C. Koutschan, A fast approach to creative telescoping, Math. Comput. Sci. 4 (2-3) (2010) 259–266.URL https://doi.org/10.1007/s11786-010-0055-0 [99] C. Koutschan, Y. Zhang, Desingularization in the q -Weyl algebra, Adv. in Appl. Math. 97 (2018)80–101.URL https://doi.org/10.1016/j.aam.2018.02.005 [100] H. Labrande, Computing Jacobi’s theta in quasi-linear time, Math. Comp. 87 (311) (2018) 1479–1508.URL https://doi.org/10.1090/mcom/3245 [101] H. Labrande, E. Thomé, Computing theta functions in quasi-linear time in genus two and above,LMS J. Comput. Math. 19 (suppl. A) (2016) 163–177.URL https://doi.org/10.1112/S1461157016000309 [102] J. Le Caine, The linear q -difference equation of the second order, Amer. J. Math. 65 (1943) 585–600.URL https://doi.org/10.2307/2371867 [103] B. Le Stum, A. Quirós, On quantum integers and rationals, in: Trends in number theory, vol. 649of Contemp. Math., Amer. Math. Soc., Providence, RI, 2015, pp. 107–130.URL https://doi.org/10.1090/conm/649/13022 [104] R. J. Lipton, Polynomials with − coefficients that are hard to evaluate, SIAM J. Comput. 7 (1)(1978) 61–69.URL https://doi.org/10.1137/0207004 [105] J.-C. Liu, Some finite generalizations of Euler’s pentagonal number theorem, Czechoslovak Math.J. 67(142) (2) (2017) 525–531.URL https://doi.org/10.21136/CMJ.2017.0063-16 [106] T. E. Mason, On Properties of the Solutions of Linear q -Difference Equations with Entire FunctionCoefficients, Amer. J. Math. 37 (4) (1915) 439–444.URL https://doi.org/10.2307/2370216 [107] J. C. P. Miller, D. J. S. Brown, An algorithm for evaluation of remote terms in a linear recurrencesequence, Comput. J. 9 (1966) 188–190.[108] S. Morier-Genoud, V. Ovsienko, On q -deformed real numbers, Exp. Math. (2019) 1–9. To appear.URL https://doi.org/10.1080/10586458.2019.1671922 [109] S. Morier-Genoud, V. Ovsienko, q -deformed rationals and q -continued fractions, Forum Math.Sigma 8 (2020) Paper No. e13, 55.URL https://doi.org/10.1017/fms.2020.9 https://doi.org/10.1090/mcom/3231 [111] C. F. Osgood, On the diophantine approximation of values of functions satisfying certain linear q -difference equations, J. Number Theory 3 (1971) 159–177.URL https://doi.org/10.1016/0022-314X(71)90033-3 [112] A. Ostrowski, On two problems in abstract algebra connected with Horner’s rule, in: Studies inmathematics and mechanics presented to Richard von Mises, Academic Press Inc., 1954, pp. 40–48.[113] I. Pak, Partition bijections, a survey, Ramanujan J. 12 (1) (2006) 5–75.URL https://doi.org/10.1007/s11139-006-9576-1 [114] V. Y. Pan, Methods of computing values of polynomials, Russian Mathematical Surveys 21 (1)(1966) 105–136.[115] M. Paterson, L. Stockmeyer, Bounds on the evaluation time for rational polynomial, in: 12thAnnual Symposium on Switching and Automata Theory (SWAT 1971), 1971, pp. 140–143.URL https://doi.org/10.1109/SWAT.1971.5 [116] M. S. Paterson, L. J. Stockmeyer, On the number of nonscalar multiplications necessary to evaluatepolynomials, SIAM J. Comput. 2 (1973) 60–66.URL https://doi.org/10.1137/0202007 [117] P. Paule, S. Radu, Rogers-Ramanujan functions, modular functions, and computer algebra, in:Advances in computer algebra, vol. 226 of Springer Proc. Math. Stat., Springer, Cham, 2018, pp.229–280.URL https://doi.org/10.1007/978-3-319-73232-9_10 [118] P. Paule, A. Riese, A Mathematica q -analogue of Zeilberger’s algorithm based on an algebraicallymotivated approach to q -hypergeometric telescoping, in: Special functions, q -series and relatedtopics, vol. 14 of Fields Inst. Commun., Amer. Math. Soc., 1997, pp. 179–210.[119] M. Petkovšek, H. S. Wilf, D. Zeilberger, A = B , A K Peters, 1996.[120] J. M. Pollard, Theorems on factorization and primality testing, Proc. Cambridge Philos. Soc. 76(1974) 521–528.URL https://doi.org/10.1017/s0305004100049252 [121] L. R. Rabiner, R. W. Schafer, C. M. Rader, The chirp z -transform algorithm and its application,Bell System Tech. J. 48 (1969) 1249–1292.URL https://doi.org/10.1002/j.1538-7305.1969.tb04268.x [122] A. Riese, qMultiSum—a package for proving q -hypergeometric multiple summation identities, J.Symbolic Comput. 35 (3) (2003) 349–376.URL https://doi.org/10.1016/S0747-7171(02)00138-4 [123] L. Rogers, S. Ramanujan, Proof of certain identities in combinatory analysis, Proc. CambridgePhilos. Soc. 19 (1919) 211–214.[124] L. J. Rogers, Second Memoir on the Expansion of certain Infinite Products, Proc. Lond. Math.Soc. 25 (1893/94) 318–343.URL https://doi.org/10.1112/plms/s1-25.1.318 [125] H. A. Rothe, Formulae de serierum reversione demonstratio universalis signis localibuscombinatorico-analyticorum vicariis exhibita, Leipzig (1793).[126] C. Sabbah, Systèmes holonomes d’équations aux q -différences, in: D -modules and microlocalgeometry (Lisbon, 1990), de Gruyter, 1993, pp. 125–147.[127] C.-P. Schnorr, Improved lower bounds on the number of multiplications / divisions which arenecessary to evaluate polynomials, Theoret. Comput. Sci. 7 (3) (1978) 251–261.URL https://doi.org/10.1016/0304-3975(78)90016-6 [128] P. Scholze, Canonical q -deformations in arithmetic geometry, Ann. Fac. Sci. Toulouse Math. (6)26 (5) (2017) 1163–1192.URL https://doi.org/10.5802/afst.1563 [129] A. Schönhage, Schnelle Multiplikation von Polynomen über Körpern der Charakteristik 2, ActaInformatica 7 (1977) 395–398.[130] A. Schönhage, V. Strassen, Schnelle Multiplikation grosser Zahlen, Computing (Arch. Elektron.Rechnen) 7 (1971) 281–292.URL https://doi.org/10.1007/bf02242355 [131] D. Shanks, A short proof of an identity of Euler, Proc. AMS 2 (1951) 747–749.URL https://doi.org/10.2307/2032076 [132] T. Sprenger, W. Koepf, Algorithmic determination of q -power series for q -holonomic functions, J.Symbolic Comput. 47 (5) (2012) 519–535. RL https://doi.org/10.1016/j.jsc.2011.12.004 [133] V. Strassen, Polynomials with rational coefficients which are hard to compute, SIAM J. Comput.3 (1974) 128–149.URL https://doi.org/10.1137/0203010 [134] V. Strassen, Einige Resultate über Berechnungskomplexität, Jber. Deutsch. Math.-Verein. 78 (1)(1976/77) 1–8.[135] H. W. L. Tanner, On the Enumeration of Groups of Totitives, Proc. Lond. Math. Soc. 27 (1895/96)329–352.URL https://doi.org/10.1112/plms/s1-27.1.329 [136] T. Tao, E. Croot, III, H. Helfgott, Deterministic methods to find primes, Math. Comp. 81 (278)(2012) 1233–1246.URL https://doi.org/10.1090/S0025-5718-2011-02542-1 [137] W. J. Trjitzinsky, Analytic theory of linear q -difference equations, Acta Math. 61 (1) (1933) 1–38.URL https://doi.org/10.1007/BF02547785 [138] H. Tsai, Weyl closure of a linear differential operator, vol. 29, 2000, pp. 747–775, Symbolic com-putation in algebra, analysis, and geometry (Berkeley, CA, 1998).URL https://doi.org/10.1006/jsco.1999.0400 [139] M. van der Put, M. F. Singer, Galois theory of linear differential equations, vol. 328 of Grundlehrender Mathematischen Wissenschaften, Springer, 2003.URL https://doi.org/10.1007/978-3-642-55750-7 [140] H. S. Wilf, D. Zeilberger, An algorithmic proof theory for hypergeometric (ordinary & q ) multi-sum/integral identities, Invent. Math. 108 (3) (1992) 575–633.URL https://doi.org/10.1007/BF02100618 [141] K.-W. Yang, On the product Q n ≥ (1 + q n x + q n x ) , J. Austral. Math. Soc. Ser. A 48 (1) (1990)148–151.[142] M. Yip, Rook placements and Jordan forms of upper-triangular nilpotent matrices, Electron. J.Combin. 25 (1) (2018) Paper 1.68, 25.[143] D. Zagier, Elliptic modular forms and their applications, in: The 1-2-3 of modular forms, Univer-sitext, Springer, 2008, pp. 1–103.URL https://doi.org/10.1007/978-3-540-74119-0_1 [144] D. Zeilberger, A holonomic systems approach to special functions identities, J. Comput. Appl.Math. 32 (3) (1990) 321–368.URL https://doi.org/10.1016/0377-0427(90)90042-Xhttps://doi.org/10.1016/0377-0427(90)90042-X