Alin Bostan

Computing the N -th Term of a q-Holonomic Sequence

Keywords: CCS CONCEPTS, Computing methodologies → Algebraic algorithms Algorithms, complexity, q-factorial, q-holonomic sequences

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 the same arithmetic complexity. We design q-analogues of 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 of any q-holonomic sequence.

Both algorithms work in arithmetic complexity quasi-linear in

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). 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 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 that for any field K, an arbitrary polynomial P ∈ K[x] of degree n can be evaluated using O ( √ n) nonscalar multiplications; 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 remains linear in n.

On the other hand, 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 repeated squaring 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 ) of univariate polynomials is called hard to compute if 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 | On the number of nonscalar multiplications necessary to evaluate polynomials[END_REF] proved the existence of polynomials in K[x] which are hard to compute. Specific families of hardto-compute polynomials were first exhibited by Strassen [START_REF] Strassen | Polynomials with rational coefficients which are hard to compute[END_REF]. 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 example is R n (x ) := n k =0 x k 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 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 / 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 theo- rem [56, §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; by [START_REF] Kirillov | On a remarkable sequence of polynomials[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 qdifferential equations, or equivalently, q-shift (or, q-difference) 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, Gauss and Heine. 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, Jackson, Carmichael, Mason, Adams, Trjitzinsky, Le Caine and Hahn, 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-difference 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-differential 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 [48, §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 q-multisums") 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-differential 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 [25, § 2.5] and [START_REF] Garoufalidis | A survey of q-holonomic functions[END_REF].

The simplest non-trivial holonomic sequence is n!, which 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!, 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 kdimensional subspaces of F n q (points on Grassmannians over F q ). There are q-analogs to (almost) everything. To select just two basic examples, the q-analog [5, Thm. 3.3] of the binomial theorem is

n k =1 (1 + q k -1 x ) = n k =0 n k q q ( k 2 ) x k (1)
and the q-version [5, 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] Pak | Partition bijections, a survey[END_REF] and other subfields of combinatorics [START_REF] Fürlinger | q-Catalan numbers[END_REF][START_REF] Kirillov | On a remarkable sequence of polynomials[END_REF]; theta functions and modular forms [START_REF] Labrande | Computing Jacobi's theta in quasi-linear time[END_REF][START_REF] Paule | Rogers-Ramanujan functions, modular functions, and computer algebra[END_REF][START_REF] Zagier | Elliptic modular forms and their applications[END_REF]; special functions [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 | Irreducibility of q-difference operators and the knot 7 4[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] Vizio | Arithmetic theory of q-difference equations: the q-analogue of Grothendieck-Katz's conjecture on p-curvatures[END_REF].

The main message of this article is that for any example of q-holonomic sequence occurring in those various fields, one can compute selected coefficients faster than by a direct algorithm.

Complexity basics. We estimate the 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 known upper bound is θ < 2.3729 [START_REF] Gall | Powers of tensors and fast matrix multiplication[END_REF]. Most arithmetic operations on univariate polynomials of degree d in K[x] can be performed in quasilinear 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 [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]. 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].

M(d ) = O (d log d log log d )

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.

De Feo's question

Here is our first example, emerging from a question asked to the author by Luca De Feo * ; 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 αq, αq 2 , . . . , αq N -1 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 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 ). (3) 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 | q-hypergeometric solutions of q-difference equations[END_REF]. Using fast polynomial multiplication, we conclude that F (α ) can be computed in arithmetic complexity quasi-linear in √ 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 babystep / 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 [62, p. 523]. 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 (1)P (Q ) • • • P (Q s-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 (α ). * Private (email) communication, 10 January, 2020. * * 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 ).
As pointed out in the remarks after the proof of [21, 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 (x ) • P t (q t x ) and

P 2t +1 (x ) = P t (x ) • P t (q t x ) • (α -q 2t 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 i j = 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 i j , 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 2M(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[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.

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 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 [56, §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 +b j 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)),
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 if q := q 2a , then the following matrix q-factorials coincide:

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

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

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

:= M ( qs-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 .

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 is q-holonomic if it satisfies a nontrivial q-recurrence, that is, a linear recurrence with coefficients polynomials in q and q n . Definition 3.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 , . . . , c r in K[x, y], with c r 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. (5)

The integer r is called the order of the q-recurrence (5). 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 . ( 6
)
They are clearly q-holonomic, and even q-hypergeometric.

The sequences

(u n ) = (q n ), (v n ) = (q n 2 ) and (w n ) = (q ( n 2
) ) 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 ) is not q-holonomic [37, Ex. 2.2(b)]. Another basic example is the q-Pochhammer symbol

(x; q) n := n-1 k =0 (1 -xq k ) (7)
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.

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, 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]. 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 (5) can be translated into a first-order recurrence on r × 1 vectors:

         u n+r . . . u n+1          =            -c r -1 c r • • • -c 1 c r -c 0 c r 1 • • • 0 0 . . . . . . . . . . . . 0 • • • 1 0            ×          u n+r -1 . . . u n          . (8)
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), (9)
where M (q n ) denotes the companion matrix from equation (8).

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 [19, §6], where the cost O (M( √ N )) is reached under some invertibility assumptions.

Theorem 3.2. 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 (α

• 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 3.2. 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). □

N -th term of a q-holonomic sequence

We give 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 [19, §6], where the improved cost O (M( √ N )) is reached under additional invertibility assumptions.

Theorem 3.4. 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] Andrews | The theory of partitions[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 ( 8), 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] Bernstein | Algorithmic number theory: lattices, number fields, curves and cryptography[END_REF]. Algorithms 3 and 4 adapt mutatis mutandis to this effect. □ Corollary 3.5. 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 truncation of degree N -1.

If α ∈ K, 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, 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.4 to v N (q) concludes the proof. □

Remark that the same result holds if e q (x ) is replaced by any power series satisfying a q-differential equation. For instance, one can evaluate fast all truncations of Heine's q-hypergeometric series

2 ϕ 1 ([a, b], [c]; q; x ) := n ≥0 (a; q) n (b; q) n (c; q) n • x n (q) n .
Remark that Theorem 3.4 can be adapted to the computation of several coefficients of a q-holonomic sequence. We omit the proof, which 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]. Theorem 3.6. Under the assumptions of Theorem 3.4, let

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

The case q is an integer: bit complexity

Until now, we only considered the arithmetic complexity model. We briefly 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, and the matrix q-factorials from §3.1 are computed by binary splitting rather than by baby-steps / giant-steps. As an illustrative example, consider the computation of the term u N (q) = N -1 n=0 q n 2 , where q is assumed to be an integer of B bits.

The integer u N (q) is bounded in absolute value by Nq N 2 , so its bitsize is of magnitude N 2 B. The "naive" algorithm consisting of computing the summands q n 2 one after the other, before summing them, has bit complexity Õ (N 3 B). This is not (quasi-)optimal with respect to the output size. Can one do better? The answer is "yes". 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 (9) as in §2.2. Now the point is that, instead of using baby-steps / giant-steps, it is a better idea to use binary splitting. The complexity of this approach becomes then quasi-optimal, that is Õ (N 2 B), which is quasi-linear in the bitsize of the output. The following general result can be proved along the same lines.

Theorem 3.7. Under the assumptions of Theorem 3.4, 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-differential 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 [26, §6], anticipated a decade earlier (without proof) by Schroeppel and Salamin in Item 178 of [7]; see [8, §12] for a good survey on binary splitting.

APPLICATIONS 4.1 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 a remarkable sequence of polynomials[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 Roger-Ramanujan identities [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 | 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 [56, §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 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-differential equations

The computation of polynomial and rational solutions of linear differential equations lies at the heart of several important algorithms, for computing hypergeometric and Liouvillian solutions, for factoring and for computing differential Galois groups [START_REF] Van Der Put | Galois theory of linear differential equations[END_REF]. Creative telescoping algorithms (of second generation) for multiple 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-differential equations: improving algorithms for polynomial and rational solutions of such equations is important in finding q-hypergeometric solutions [START_REF] Abramov | q-hypergeometric solutions of q-difference equations[END_REF], in computing q-differential Galois groups [START_REF] Hendriks | An algorithm for computing a standard form for second-order linear q-difference equations[END_REF], and in 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. The same improvements can be transferred to q-differential 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 * .

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 [61, §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-differential case, Koornwinder's q-Gosper algorithm [49, §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

A preliminary implementation in Magma of Algorithms 1 and 2 in §2.1 delivers 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. In our experiments, K is the finite field F p with p = 2 30 + 3 elements. Timings are given 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: timings are 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.

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-differential equations.

Here are some questions that we plan to investigate in the future. 1. (Computing curvatures of q-differential equations) In the differential case, p-curvatures can be computed 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 pcurvature[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]. What about the q-differential analogue? One strong motivation comes from the fact that the q-analogue [10] of Grothendieck's conjecture (relating equations over Q with degree N Naive algorithm Algorithm 1 Algorithm 2 2 16 0.04 0.03 0.00 2 18 0.18 0.03 0.01 2 20 0.72 0.06 0.01 2 22 2.97 0.14 0.02 2 24 11.79 0.32 0.04 2 26 47.16 0.73 0.08 2 28 188.56 1.68 0.15 2 30 755.65 3.84 0.31 2 32 3028.25 8.65 0.64 2 34 19.65 1.41 2 36 44.42 2.96 2 38 101.27 6.36 2 40 228.58 14.99 2 42 515.03 29.76 2 44 1168.51 61.69 2 46 2550.28 137.30 2 48 297.60 2 50 731.63 2 52 1395.33 2 54 3355.39

Table 1 Comparative timings (in seconds) for the computation of N -1 i =0 (α -q i ) ∈ F p , with p = 2 30 + 3 and (α, q) randomly chosen in F p × F p . All algorithms were executed on the same machine, running Magma v. 2.24. For each target degree N , each execution was limited to one hour. Naive algorithm could reach degree N = 2 32 , Algorithm 1 degree N = 2 46 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]. This could be used to improve the computation of rational solutions. 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 2k x 2 +aq k x +b). This might have applications related to Question 1, or to counting points on q-deformations [67]. 3. (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. How fast can one compute (truncations / evaluations of) quantized versions of numbers like e or π ? 4. (Evaluating more polynomials) Is it possible to evaluate fast polynomials of the form N n=0 x n 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 in isogeny-based cryptography.

  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[START_REF] Abramov | On polynomial solutions of linear operator equations[END_REF]. By proceeding as in Algorithm 2 in §2.1, the complexity of Algorithm 3 already is quasi-linear in √ N . However, its dependence in a, b 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 M (x ) := prx + 1 -prx 1 0 , with r := q a+b ,

- 1 )Corollary 3 . 3 .

 133 in O (M( √ N )) operations in K. The cost of computing r N and α N is O (log N ), and thus it is negligible. □ Under the assumptions of Theorem 3.2 and for any n ∈ N, one can compute in O (M( √ N )) operations in K:

  , 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. (Computing points on q-curves) Counting efficiently points on (hyper-)elliptic curves leads to questions like: for a, b ∈ Z, compute the coeff. of x p-1 2 in G p (x ) := (x 2 + ax + b) p-1 2

* As usual, the notation Õ (•) is used to hide polylogarithmic factors in the argument.

* A technical subtlety is that, as pointed out in [1, §4.3], it is not obvious in the qdifferential 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[22, §3.2].

Acknowledgements. I thank Luca De Feo for his initial question, who motivated this work, and for the very interesting subsequent discussions. My friendly thanks go to Lucia Di Vizio, Kilian Raschel and Sergey Yurkevich for their careful reading of the manuscript. I am indebted to the three referees for many helpful remarks. This work was supported in part by DeRerumNatura ANR-19-CE40-0018.


Djamel Ait-Akli

Abdelkader Merakeb

A Trace Result for Sobolev Functions on Planar Jordan Domains

Keywords: Sobolev spaces, Jordan domains, trace inequality, density of smooth functions, Lipschitz regularity, Hausdorff measure AMS subject classification: 46E35, 26A16

After recalling some facts from Hausdorff measure theory, we define, using the notion of s-length, the space L p (∂Ω, d σ s ) on the boundary of a Jordan domain Ω. Then, with the aid of some auxiliary considerations, we establish a trace result for Sobolev's space functions W 1,p (Ω), 1 < p < ∞ in case Ω ⊂ R 2 is an arbitrary bounded Jordan domain and ∂Ω has a finite s-measure with s > 1. This result generalizes the classical case for Lipschitz regular domains.

Introduction

Let Ω ⊂ R2 be a bounded Jordan domain i.e. its boundary, denoted ∂Ω, is homeomorphic to the boundary, ∂B 1 , of the unit ball, B 1 , in R 2 . We denote diam Ω the diameter of Ω. Let us note W k,p (Ω) the usual Sobolev spaces defined by

W k,p (Ω) := {v ∈ L p (Ω) such that ∂ α v ∈ L p (Ω) ∀α, multi-index, |α| ≤ k}.
The trace operator, applied to functions of the space W1,p (Ω) defined on a domain with Lipschitz boundary, is a well-established notion in Sobolev's space theory. We propose to show that a trace operator can also be defined when the domain is characterized by a lesser regularity, more precisely, when it is only Jordan. An important fact about the class of Jordan domains is that they generally do not admit the property of (1, p)-extension, see [START_REF] Maz'ya | Extension of functions from Sobolev spaces, English translation[END_REF]. The property of (1, p)-extension domain Ω means that an extension operator can be defined from Sobolev space W 1,p (Ω) to W 1,p (R 2 ), we can refer to [START_REF] Lewis | Approximation of Sobolev functions in Jordan domains[END_REF] for a detailed discussion of this class of domains. Lewis has showed, see [START_REF] Lewis | Approximation of Sobolev functions in Jordan domains[END_REF], that C ∞ (Ω) is dense in W 1,p (Ω) when Ω ⊂ R 2 is an arbitrary planar Jordan domain. We establish, using, among other things, this fact, a trace theorem for the functions in W 1,p (Ω) on the boundary Γ := ∂Ω of the Jordan domain Ω.

Let us Recall the continuity property of the trace operator on the boundary of a Lipschitz (or weakly Lipschitz) regular domain D ⊂ R 2 defined for the Sobolev space W 1,p (D), see ( [START_REF] Gagliardo | Caratterizzazioni delle tracce sulla frontiera relative ad alcune classi di funzioni in n variabili[END_REF]): Proposition 1. Let D ⊂ R 2 be a domain with Lipschitz boundary and 1 < p < ∞. The space W 1,p (D) is continuously embedded into L p (∂Ω) i.e. there exists a constant C D > 0 such that for all u ∈ W 1,p (D) we have

||u|| p,∂Ω ≤ C D ||u|| W 1,p (D) . (1)
It should be noted that proposition [START_REF] Lewis | Approximation of Sobolev functions in Jordan domains[END_REF] does not apply to our Jordan domain Ω given the assumption of Lipschitz regularity of the domain D. Let us also recall the density result proved in ( [START_REF] Lewis | Approximation of Sobolev functions in Jordan domains[END_REF], p.261):

Theorem 2. If Ω ⊂ R 2 is a domain whose boundary is a Jordan curve then C ∞ (Ω) is dense in W 1,p (Ω) for 1 < p < ∞.
We present hereafter some facts from Hausdorff measure theory, this later constitutes the point of view under which the corpus of the paper is to be considered.

Hausdorff measure

We recall the definition of the Hausdorff measure of a set. We can refer to ( [START_REF] Yeh | Real Analysis, Theory of measure and integration[END_REF], chap.7) for the results and facts stated in this sub-section. Let A ⊂ R 2 . We call a δ-cover of the subset A, the collection

R δ (A) := {E i ⊂ R 2 | diam E i | ≤ δ A ⊂ ∪ i E i },
and we define the exterior, auxiliary, measure

H s δ (A) := inf{ i∈N | diam E i | s such that E i ∈ R δ }. (2)
The Hausdorff s-measure of A is then defined by

H s (A) := lim s→0 H s δ (A). (3)
It's helpful to emphasize, see [START_REF] Yeh | Real Analysis, Theory of measure and integration[END_REF], that the measure

H s δ is decreasing in δ i.e. if A is a set then 0 < δ 1 < δ 2 implies H s δ 2 (A) < H s δ 1 (A). (4)
On the other hand, the Hausdorff measure is an increasing set function in the sense that if A and B are two sets then

A ⊂ B implies H s (A) ≤ H s (B). (5)
Lastly, the Hausdorff dimension, dim H (A), of a set A is defined by

dim H (A) := inf s>0 {H s (A) = 0} = sup s>0 {H s (A) = +∞}. ( 6
)
We recall the definition of Hausdorff's distance. Let (G, d) be a metric space and G F be the collection of non-empty bounded closed subsets of G. The Hausdorff distance on

G F is the application dist H , of G F ×G F in R + , defined by (X, Y ) → dist H (X, Y ) := max{sup x∈X inf y∈Y d(x, y), sup y∈Y inf x∈X d(x, y)}. ( 7
)
Here and in the sequel, the distance (x, y) → d (x, y) is understood in the sense of the euclidean distance between the points x and y.

Remark 3. 1)-Given a Jordan domain Ω ⊂ R 2 , we can construct a nested sequence of Lipschitz domains (Ω n ) n∈N that approximates Ω in the sense of Hausdorff 's distance [START_REF] Taylor | Measure Theory and Integration[END_REF] with Ω n ⊂⊂ Ω n+1 ⊂ Ω. Therefore, the Jordan curve Γ = ∂Ω can be approximated by a sequence of Lipschitz curves Γ n .

In addition, we choose a precise parameterization γ m for each of the curves Γ m ; these functions are defined by

γ m : [a, b] → Γ m t → γ m (t) (8)
such that γ m is defined as follows: if γ(t) = (γ 1 (t), γ 2 (t)) is the point of Γ corresponding to t then the point γ m (t) is the unique point of

Γ m that satisfies |γ m (t) -γ(t)| 2 = inf x∈Γm d (x, γ(t)) .
2)-We should also recall that the n-dimensional Hausdorff measure for n ∈ N * coincides with Lebesgue's n-dimensional measure up to a multiplicative constant C λ n .

We assume, in the sequel, that the domain Ω is Jordan with boundary a Jordan curve Γ = ∂Ω ⊂ R 2 which is Hausdorff s-dimensional with s > 1 such that ∂Ω = Γ has finite s-dimensional measure. The main result of this paper is stated in the following theorem:

Theorem 4.
Let Ω be Jordan domain and p such as 1 < p < ∞. We assume that the boundary ∂Ω has finite s-measure, i.e. H s (∂Ω) < ∞. Then, we can define a trace operator

T : W 1,p (Ω) → L p (∂Ω) u → T u
which extends the usual restriction operator to ∂Ω for continuous functions. Furthermore, T is continuous i.e. there is a constant c t > 0 independent of u such that for every u ∈ W 1,p (Ω)

||T u|| p,∂Ω ≤ C t ||u|| W 1,p (Ω) , (9)
where we have noted || || p,∂Ω the L p space according to a measure that we will define below and that is related to the Hausdorff measure on the boundary.

Stieltjes integral and s-length

Let γ be a continuous, C 0 ([a, b]), parametric representation of Γ. Since we are dealing with a curve Γ of Hausdorff dimension s > 1 then the usual arc length, i.e. with respect to Lebesgue's measure or, similarly, to Hausdorff 1-measure, as well as the curvilinear abscissa are obsolete in our situation since they may equal +∞. This suggests adapting the definition of the curvilinear abscissa to our situation. The curvilinear s-abscissa is defined by

t → s γ (t) := H s (γ([a, t])); (10)
i.e. the Hausdorff s-dimensional measure of the "piece" of the curve that is the image of [a, t] by γ. It is also the s-length of the arc γ restricted to the interval [a, t], see Definition 5 below. We define the curvilinear integral of the continuous function u : Γ → R as the Stieltjes integral of u•γ with respect to the curvilinear s-abscissa defined by (10). The notion of curvilinear s-abscissa is the key to understanding this integral. Indeed, the Stieltjes integral in question is given by:

Γ u d s = [a,b] u • γ(t) d s γ (t)
where the integral can be, at first, intuitively understood in the sense of Riemann; i.e. it is the limit, when the maximal step σ := max i [t k , t k+1 ] of the partition of [a, b] tends towards zero, of the Riemann sums:

lim σ→0 u(γ(t k )) (s γ (t k ) -s γ (t k-1 )) (11)
where the partition is denoted :

a = t 0 < t 1 < . . . < t n = b, t k ∈ [t k-1 , t k ].
This definition does not depend on the parameterization chosen for Γ. Note that the limit (11), and thus the integral, exists as soon as u is bounded; this in addition to the fact that the curvilinear s-abscissa is of bounded variation; indeed, since ]0,

t k ] ⊂]0, t k+1 ] implies γ(]0, t k ]) ⊂ γ(]0, t k+1 ]
) and that H s is increasing in the sense of (11), then the curvilinear s-abscissa function is increasing and bounded. This also implies that its derivative,

s γ (t
interval ]u, v] ⊂ [a, b] w(]u, v]) = s γ (v) -s γ (u) = H s (γ(]0, v])) -H s (γ(]0, u])),
by the Caratheodory extension theorem, there is a unique Borel measure µ s on [a, b] equal to w on any interval. The measure µ s emerges from the exterior measure defined by

µ s (E) = inf i w(I i ) : E ⊂ i I i (12)
the infimum being taken over the coverings of E by countable semi-open intervals. This measure is called the Lebesgue-Stieltjes measure associated with s γ . This measure generalizes the Lebesgue measure in the extent that it coincides with the latter in the particular case when w(]u, v]) = v -u. More rigorously, we say that a bounded or continuous function f is integrable with respect to the Stieltjes measure or Stieltjes integrable if f is integrable in the Lebesgue-Stieltjes sense i.e. if we have b a |f (t)| ds γ (t) < ∞ which, in turn, is defined as Lebesgue's classic integral of |f | with respect to the measure µ s . In particular the Stieltjes integral of f coincides with the integral of the function f s γ in Lebesgue's sense with respect to the Lebesgue measure. In our situation s γ exists almost everywhere. On the other hand, we have not explicitly expressed the dependency of the measure µ s and later that of µ s δ on the parameterization γ but it is to be understood implicitly, in the sequel, that the measures µ s and µ s δ are defined according to the context relatively to Γ or Γ n , the same remark holds for the set function that we denoted w instead of w γ .

Definition 5. Let Γ be a Jordan curve in R 2 and γ : [a, b] → R 2 a continuous, C 0 ([a, b]
), parametric representation of Γ. We assume that Γ is an s-dimensional Hausdorff set. We define the s-length of Γ to be L s (Γ) := sup P,partition

V s I (γ, P ) (13)
with

V s I (γ, P ) = np-1 k=0 d s (γ(t k+1 ), γ(t k )) ( 14
)
where the supremum browses all the partitions

P := {p = {t 0 , t 1 , ..., t np }|p is a partition of the interval [a, b]}.
If L s (Γ) < ∞ then we say that Γ is s-rectifiable and that t → γ(t) has s-bounded variation.

Lemma 6. Let Γ be a Jordan curve in R 2 such as in definition ( 5), parametrized by

γ ∈ C 0 ([a, b]) then H s (Γ) = L s (Γ).
That is, the s-length of Γ equals its s-Hausdorff measure.

Proof. The proof is done in two steps 1) We show that H s (Γ) ≤ L s (Γ). We construct inductively a partition of [a, b] as follows:

t 0 := a; t k+1 := inf{x ∈ [t k , b]; d s (γ(t k ) -γ(x)) = } ∪ {b},
so we get a partition of [a, b] defined by P := {t 0 , t 1 , ..., t np } such as

d s (γ(t k ), γ(t k+1 )) = diam s (γ([t k , t k+1 ])) ∀k ∈ {0, ..., n p -2} d s (γ(t np ), γ(t np-1 )) = diam s (γ([t np , t np-1 ])) ≤ .
We have then:

H s (γ) ≤ np-2 k=0 diam s (γ([t k , t k+1] )) + diam s (γ([t np , t np-1] )) ≤ np-2 k=0 d s (γ(t k ), γ(t k+1 )) + ≤ L s (Γ) + , letting → 0 we find H s (Γ) ≤ L s (Γ).
2) We show that L s (Γ) ≤ H s (Γ). It is easy to see that:

d (γ(t k+1 ), γ(t k )) ≤ diam(γ([t k-1 , t k ])) ≤ H 1 (γ([t k-1 , t k ])), thus d s (γ(t k+1 ), γ(t k )) ≤ H 1 (γ([t k-1 , t k ])) s .
Furthermore, for all > 0, there exists 0 < δ < 1 such that

d s (γ(t k+1 ), γ(t k )) ≤ H 1 δ (γ([t k-1 , t k ])) s + .
According to Definition 2 of the exterior measure

H 1 δ , if (A i ) i is such that A i ∈ R δ is a δ -cover of γ([t k-1 , t k ])) then we have d s (γ(t k+1 ), γ(t k )) ≤ [ i diam(A i )] s + .
Since diam(A i ) < 1 for all i then we get

d s (γ(t k+1 ), γ(t k )) ≤ i diam s (A i ) + , for all A i ∈ R δ such that diam(A i ) ≤ δ . We then deduce that d s (γ(t k+1 ), γ(t k )) ≤ inf R δ i diam s (A i ) + .
Letting → 0, we deduce, using ( 3) and (4) that

d s (γ(t k+1 ), γ(t k )) ≤ H s (γ([t k+1 ) -γ(t k )]).
Since H s is a sigma-additive measure, then

L s (Γ) = d s (γ(t k+1 ), γ(t k )) ≤ H s (Γ).
The lemma (6) makes it clear that the integral of the constant function 1, which by definition of the integral is equal to H s (Γ), is nothing else than the s-length of Γ defined by ( 13)-( 14). This fact is present especially in the case of the dimension s = 1 where the 1-length is nothing other than the usual length of the curve. One can also define (or construct) the curvilinear integral with respect to the Hausdorff s-dimensional measure, with which the curvilinear s-abscissa is expressed on Γ. This will be the subject of the next sub-section. We emphasize that an almost similar work has been presented in the context of Lebesgue-Hausdorff line integral, see ( [START_REF] Lazwardi | Lebesgue-Hausdorff Line Integral of Hausdorff Measurable Multivariable Function over Simple Curve on[END_REF]).

Construction of the Lebesgue-Hausdorff-Stieltjes (LHS) line integral on the Jordan curve

Let γ : [a, b] → Γ, γ ∈ C 0 ([a, b]
), be parametrization of the Jordan curve Γ. Let ψ be a simple defined on [a, b] i.e. there exists a sequence of real numbers (α i ) i and a collection of sets

(E i ) i such that [a, b] = ∪ i E i and ψ(t) = i α i χ(E i ).
Then, the integral of ψ is defined by

(LHS) [a,b] ψ d µ s = i α i µ s (E i ),
where the measure µ s is defined by (12). The lower integral of the continuous function u • γ defined on [a, b] is given by

(LHS) Γ u d σ s = [a,b] u•γ d µ s = sup{ [a,b] ζ d µ s , ζ ≤ u•γ , ζ simple function on [a, b]}.
The upper integral of the function u • γ is defined by

(LHS) Γ u d σ s = [a,b] u•γ d µ s = inf{ [a,b] ξ d µ s , u•γ ≤ ξ , ξ simple function on [a, b]}. Now, the function u • γ is said to be s-Hausdorff-Lebesgue-Stieltjes inte- grable if (LHS) [a,b] u • γ d µ s = (LHS) [a,b] u • γ d µ s ,
this number is the line integral of u on Γ and is denoted

(LHS) Γ u d σ s .
If the function |u • γ| is integrable then we note u ∈ L 1 (Γ, d σ s ). We define also the space L p (Γ, d σ s ) as the set of measurable functions u such that |u| p ∈ L 1 (Γ, d σ s ). |u • γ| p dµ s δ can be defined in exactly a similar way as done above for the measure µ s .

2)-On the other hand, using the definition of the integral with respect to the measure µ s δ , one can easily see that

lim δ→0 Γ |u| d σ s δ = lim δ→0 [a,b] |u•γ(t)| d µ s δ (t) = [a,b] |u•γ(t)| d µ s (t) = Γ |u| d σ s .
(15) 3)-Assume that u is Lipschitz regular on Ω then consider the functions u • γ and u • γ m as defined on the measure space ]a, b[,B(]a,b[),µ 1 δ ; they are clearly µ 1 δ -integrable. Thus, the definition of the integral with respect to the measure µ 1 δ and Definition 8 of the functions γ m , see Remark 3, yield, using the dominated convergence theorem applied on the sequence of functions (u • γ m ) m , the following limit:

| [a,b] |u • γ m (t)| p d µ 1 δ (t) - [a,b] |u • γ(t)| p d µ 1 δ (t)| → 0 as m → ∞. (16)
4)-We can choose m such that Γ m and Γ are arbitrarily close in the sense of the Hausdorff distance [START_REF] Taylor | Measure Theory and Integration[END_REF] i.e. such that to make

dist H (Γ, Γ m )
small enough; as a consequence we have: given 0 < δ < 1

|H s δ (Γ m ) -H s δ (Γ)| → 0 as m → ∞. (17)
This implies in particular, knowing that H 1 δ (Γ) < ∞, that there exists m δ such that for all m > m δ we have

H 1 δ (Γ m ) ≤ 2H 1 δ (Γ). ( 18
)
2 Proof of the main result

We now turn to the proof of Theorem 4.

Proof. Let p > 1 and fix δ, 0 < δ < 1. According to -1) of Remark 3, there exists a sequence of Lipschitz curves, denoted Γ n := ∂Ω n , that approximate the Jordan curve Γ := ∂Ω. Define F := L p (∪ n Γ n , d σ 1 δ ) to be the vector space whose elements are functions u defined on ∪ n Γ n such that if the parametrization functions, γ n , are as defined by (8) then, |u| is L p -integrable on the countable union, ∪ n Γ n , of the subsets (Γ n ) n with respect to the measure σ 1 δ . It is clearly a normed vector space when equipped with the norm

||v|| F := ∪nΓn |v| p d σ 1 δ 1 p = n [a,b] |v • γ n | p d µ 1 δ 1 p
, where the measure µ s δ is defined by (12). Pose E := W 1,p (Ω) ∩ C ∞ (Ω) and

F := L p (∪ n Γ n , d σ 1 δ ).
Define the sequence of linear maps (L m ) m∈N by

L m : (E, || || W 1,p (Ω) ) → (F, || || L p (∪nΓn,d σ 1 δ ) ) u → L m u with (L m u)(x) = u(x) , x ∈ Γ m 0 , x ∈ ∪ n Γ n -Γ m .
On one hand, (L m ) m is a sequence of continuous linear maps. Indeed; using the properties of the auxiliary measure H 1 δ , namely (3), ( 4) and (12), we have

||L m (u)|| p F = Γm |u| p d σ 1 δ = [a,b] |u • γ m | p d µ 1 δ ≤ [a,b] |u • γ m | p d µ 1 , (19)
using 2) of remark 3 we write

[a,b] |u • γ m | p d µ 1 = C λ 1 [a,b] |u • γ m | p d λ 1 , (20)
where λ 1 is the Lebesgue 1-dimensional measure. According to estimate (1) of the trace inequality on the Lipschitz curve Γ m = ∂Ω m we have: for all m > 0 there exists C m > 0 such that for all u ∈ W 1,p (Ω) ∩ C ∞ (Ω)

||L m (u)|| F ≤ C λ 1 [a,b] |u • γ m | p d λ 1 1 p (21) ≤ (C λ 1 ) 1 p C m ||u|| W 1,p (Ωm) ≤ (C λ 1 ) 1 p C m ||u|| W 1,p (Ω) ;
the constants C m , which are respectively the norms of the trace operators relatively to the Lipschitz domains Ω m , are finite for every m. So, using the definition of the norm of the linear maps L m , ∀m, ∃C m ≤ (C λ 1 )

1 p C m such that C m = ||L m || L(E,F ) = sup u∈E ||L m (u)|| F ||u|| W 1,p (Ω) ,
thus, for all m, there exists C m > 0 such that for all u ∈ E we have

||L m (u)|| F = [a,b] |u • γ m (t)| p d µ 1 δ (t) 1 p ≤ C m ||u|| W 1,p (Ω) . (22)
On the other hand, let u ∈ W 1,p (Ω) ∩ C ∞ (Ω) be an arbitrary fixed element. Then (18) yields us, for the fixed δ, with the following sup

m>m δ ||L m (u)|| F ≤ |u| ∞,Γm sup m>m δ [H 1 δ (Γ m )] 1 p ≤ 2|u| ∞,Ω [H 1 δ (Γ)] 1 p ≤ N (δ) < ∞. (23)
Actually, the operator L m depends also on the fixed δ, but We have not made this dependence explicit for clarity of writing and one should well notice that N (δ) → ∞ when δ → 0. Estimates ( 22) and (23) ensure, using the Banach Steinhaus theorem, the existence of a constant M (δ) > 0 which is independent of m but depending on δ such that sup m>m δ ||L m || ≤ M (δ). This implies that C m < M (δ) for all m > m δ . Therefore, there exists

C (δ) > 0 such that C m → C (δ) when m → ∞. Let u ∈ W 1,p (Ω) ∩ C ∞ (Ω).
Letting m → ∞ in (22), and using ( 16), yields us

[a,b] |u • γ(t)| p d µ 1 δ (t) 1 p ≤ C (δ)||u|| W 1,p (Ω) . (24)
The facts s > 1 and δ < 1 in addition to formulas (2) and (3) defining the measures give us

Γ |u| p d σ s δ 1 p = [a,b] |u • γ(t)| p d µ s δ (t) 1 p ≤ [a,b] |u • γ(t)| p d µ 1 δ (t) 1 p
.

(25) Combining estimates (24) and (25) yields using the Holder inequality

∀u ∈ W 1,p (Ω) ∩ C ∞ (Ω) , Γ |u| p d σ s δ 1 p ≤ C (δ)||u|| W 1,p (Ω) , (26)
the constants C (δ) are finite for every δ < 1. Define the sequence of linear maps (L δ ) δ<1 by

L δ : W 1,p (Ω) ∩ C ∞ (Ω), || || W 1,p (Ω) → W 1,p (Ω) ∩ C ∞ (Ω), || || L p (Γ,d σ s δ ) u → L δ u = u.
Estimate (26) shows clearly that the linear operator L δ is continuous for every 0 < δ < 1. Additionally L δ is obviously one-to-one. We deduce using the Banach isomorphism theorem that the inverse operator L -1 δ is also continuous i.e. there exists C -1 (δ) > 0 such that for all u ∈ W 1,p (Ω)∩C ∞ (Ω) we have

||u|| W 1,p (Ω) ≤ C -1 (δ) Γ |u| p d σ s δ 1 p , (27)
using the estimate (3) and (4) defining the auxiliary measure, we deduce immediately that

∀u ∈ W 1,p (Ω) ∩ C ∞ (Ω) , ||u|| W 1,p (Ω) ≤ C -1 (δ 0 ) Γ |u| p d σ s 1 p (28)
for some 0 < δ 0 < 1. Using again the Banach isomorphism theorem, similarly as done above, yields us, by using (28), the existence of a constant

C t > 0 such that ∀u ∈ W 1,p (Ω) ∩ C ∞ (Ω) , Γ |u| p d σ s 1 p ≤ C t ||u|| W 1,p (Ω) . (29)
According to Theorem 2, there exists (u l ) l∈N with u l ∈ C ∞ (Ω) and such that

||u l -u|| W 1,p (Ω) → 0 as l → +∞. ( 30
consequently (u l ) l is a Cauchy sequence in the complete space L p (Γ, d σ s ). Using (30), we pass into the limit in (31) to conclude that there exists u * ∈ L p (Γ, d σ s ) such that

Γ |u * | p d σ s 1 p ≤ C t ||u|| W 1,p (Ω) . (32)
We define the operator T : W for all u ∈ W 1,p (Ω); this does make T a continuous trace operator on the Jordan curve ∂Ω.

One can eventually define the space W 1-1 p ,p (∂Ω) as the range of the operator T . Thus, once its range is restricted to W 1-1 p ,p (∂Ω), the operator T is surjective.

One should notice that the continuity constant of the trace operator, C t , is, somewhat, directly related to the norms of the trace operator respectively defined on the boundary of Ω m up to the constant C λ 1 .

Conclusion

In addition to being interesting in itself, this result, like any other trace theorem for Sobolev space functions, is strongly expected to have applications in the field of mathematical analysis of partial differential equations given its wide and extensive use in the analysis literature. In particular, this type of result has proved to be very useful in demonstrating the existence and uniqueness results for PDEs.

On the other hand, it seems somewhat difficult to further generalize this result to the case p < 1, at least in the context of the above presented proof. This is due, in part, to the fact that the essential ingredient of the proof, the density result which is the subject of the theorem (2), may not be valid in this general case. Also, L p is no longer a Banach space for p < 1.

Remark 7 .

 7 1)-The line integral with respect to the measure µ s δ Γ |u| p d σ s δ = [a,b]

  ), exists almost everywhere in [a, b]. The curvilinear s-abscissa, s γ , is a set function which in turn defines a measure. This measure is called the Lebesgue-Stieltjes measure associated with s γ . To see it, we define the set function w on the semi-open generic

L2CSP, UMMTO. 15000, Tizi-Ouzou, Algeria. (djamel.aitakli@ummto.dz)

L2CSP, UMMTO. 15000, Tizi-Ouzou, Algeria.(kader.merakeb@ummto.dz)


Kevin Loquin
email: kevin.loquin@gmail.com

Isabelle Bloch
email: isabelle.bloch@telecom-paristech.frkn@cicoph.org

Kiyoko Nakashima

Florence Rossant
email: florence.rossant@isep.frmp@cicoph.org

Pierre-Yves Boëlle

Michel Paques

Pierre-Yves Boelle

Automatic Photoreceptor Detection in In-Vivo Adaptive Optics Retinal Images: Statistical Validation

Keywords: Adaptive Optics, Photoreceptor detection, in vivo diagnosis, retina imaging

Introduction

Adaptive Optics (AO) is a new ophtalmic imaging system inherited from spectacular technological advances in astronomy [START_REF] Roddier | Adaptive Optics in Astronomy[END_REF]. An AO system is made up of a deformable mirror which compensates (it can be viewed as some kind of physical deconvolution) for the image defects due to atmospheric distortions. The imaging of the retina with an AO system [START_REF] Roorda | Adaptive optics scanning laser ophthalmoscopy[END_REF] is particularly interesting because the light distortions due to the iris, pupil and lens system of the eye are similar to the distortions due to atmosphere.

The AO imaging modality technically outperforms any other up-to-date retinal imaging modality. It allows the imaging of the living human retina at the microscopic scale. For instance, the rtx1 AO system has a resolution of 0.8µm per pixel, while usual fundus cameras reach resolutions from 6 to 10 µm per pixel. This feature is illustrated in Figure 1. This new imaging modality is of prime importance for a clinical use in routine. At a microscopic level, the retina is a stack of several neuronal layers whose aim is to transform the incident light in an electrical signal that will be transmitted to the brain [START_REF] Curcio | Human photoreceptor topography[END_REF]. Within this process, the photosensitive neurons are the photoreceptor cells (rods and cones). Many inherited and acquired diseases or disorders provoke degeneration of the retina [START_REF] Ryan | Retina (3 Volume Set)[END_REF][START_REF] De | Age-related macular degeneration[END_REF], a progressive loss of the cones cells. It leads to severe visual handicap. In some diseases such as retinitis pigmentosa, visual loss can occur very early in life. In [START_REF] Wolfing | High-resolution retinal imaging of cone-rod dystrophy[END_REF], the authors illustrate the diagnosis power of retinal AO images on degenerative cones distrophy by comparing visual diagnosis procedures from these AO images (like cone counting, manual delineation of degenerated area,...) with well managed usual diagnosis procedures performed with other modalities.

The algorithm that we present in this article is an extension of an algorithm already proposed by the authors [START_REF] Loquin | Photoreceptor detection in in-vivo adpative optics images of the retina: towards a simple interactive tool for the physician[END_REF]. This algorithm is based on a recursive construction of thresholded connected components whose seeds are the regional maxima of the image. The stopping rules of this recursive construction are deduced from geometrical arguments: preference is given to brighter spots for closed merged spots which are thus considered as noisy spots. In this present article we propose to use an additional preliminary deconvolution step. The deconvolution of AO images has been proved to be particularly useful [START_REF] Blanco | Marginal blind deconvolution of adaptive optics retinal images[END_REF]: it considerably facilitates the distinction between background and information (cones) on retinal AO images.

Rare other automatic photoreceptors detection algorithms have already been proposed in the literature [START_REF] Li | Automated identification of cone photoreceptors in adaptive optics retinal images[END_REF][START_REF] Xue | Photoreceptor counting and montaging of enface retinal images from an adaptive optics fundus camera[END_REF]. Both are based on the detection of regional maxima; in [START_REF] Li | Automated identification of cone photoreceptors in adaptive optics retinal images[END_REF], the image is pre-filtered and post morphologically processed and in [START_REF] Xue | Photoreceptor counting and montaging of enface retinal images from an adaptive optics fundus camera[END_REF], the regional maxima are processed in a decreasing intensity order to discriminate the dimest ones. In this present article, we propose a validation of our algorithm by comparing it with the ones presented in [START_REF] Li | Automated identification of cone photoreceptors in adaptive optics retinal images[END_REF][START_REF] Xue | Photoreceptor counting and montaging of enface retinal images from an adaptive optics fundus camera[END_REF].

In Section 2, we present our photoreceptor detection algorithm. Then, Section 3 details the performed experiments to validate our tool, before concluding.

Photoreceptor detection algorithm

Previous version [5]

In the used AO images, the photoreceptors that we aim at detecting are bright spots, i.e. pixels with high grey levels. Hence a first natural step of our algorithm is a simple regional maxima detection. Afterwards, a maximum area of a regional maximum is defined up to a tolerance parameter, denoted by T . This tolerance parameter finds its roots in the component tree decomposition of an image [START_REF] Couprie | Quasi-linear algorithms for the topological watershed[END_REF]. In this scope, an image is seen as a decomposition of α-level connected components. The level range corresponding to the greyscale range of the image, an image is decomposed in a stack of connected components. The tolerance parameter T that we are using in our algorithm tunes the height of the top part of this stack that we consider as a maximum area.

Note that the proposed algorithm is a truncated adaptation of the watershed by immersion algorithm [START_REF] Vincent | Watersheds in digital spaces: an efficient algorithm based on immersion simulations[END_REF] when applied to the inverted image. In the inverted image, the origins of the basins used in the watershed by immersion are the regional maxima of the original image. In our case, we are not interested in reaching a complete separation of the basins (the maxima) but only in elicitating these basins until a given height, i.e. the maximum areas up to a tolerance parameter.

What's new?

Even for adaptive optics which is physically corrected by means of deformable mirrors, an acquisition necessarily transforms the real information into an altered one: a measurement. The relationship between the physical reality and its measurement is supposed to be linear, which leads to model any acquisition by a convolution operator. Deconvolution of a signal is, in some sense, reversing the acquisition process to reach as closely as possible the physical reality. With AO retinal images, the results of the deconvolution are particularly interesting. In [START_REF] Blanco | Marginal blind deconvolution of adaptive optics retinal images[END_REF], the authors propose to learn the parameters of the deconvolution kernel, i.e. the point spread function (PSF) of the sensor model, from the entire image: this is called blind deconvolution. Their results are impressive but their algorithm is too much time consuming in a clinical context. In our case, we propose a simpler approach, less computationally demanding: we use a simple Gaussian convolution model of the PSF with a unique parameter: σ, its standard deviation.

Figure 2 illustrates the interest of this technique. The right column shows (for illustration) the obtained deconvoluted images from the left column images. On the left column we can also observe the results of our method for the following 1a 1b 2a 2b The results seem to be as good as with the blind deconvolution method. However, we add a new parameter to our method, which makes the physician task more complex. In this illustration, the parameters are chosen to reach the best possible results in a subjective judgment.

Algorithm implementation: an interactive tool for the physicians

This algorithm is very rapid: after a complete cover of the image for finding the regional maxima, the recursive coverage is only applied to these extracted spots which highly diminishes the computational cost. When java-implemented as an imageJ plugin (which is not the most rapid computational framework) the algorithm result is quite instantly displayable. This remark is very important because it allows us to provide an interactive tool to the final user, the physician. This non algorithmic part of our work is greatly appreciated by the physicians who can thus easily interact with the algorithm and adapt its parameters to his visual expertise in an intuitive way. This feature highly contributes to the fact that, even if still in development, this tool is already in use in a clinical context. Another way to interact on the photoreceptors detection is to incorporate in our imageJ plugin the possibility to add non automatically detected cones and to remove wrongly detected cones manually.

Validation

Material and methods

The AO images used along this validation procedure are obtained with two prototypic systems: INOVEO and rtx1 manufactured by Imagine Eyes1 , Orsay, France. Those systems are currently operational in clinical and research settings at the Clinical Investigation Center 503 of the XV -XX hospital. According to the used system (INOVEO and rtx1), in each image, for an axial length of 24.7mm, one pixel corresponds respectively to a 1.2µm × 1.2µm and 0.8µm × 0.8µm spot on the retina. The gold standard dataset is constructed from 52 patch images of 13 bigger AO images. Thus, for each of the 13 images, 4 patch images were extracted. The patch image sizes were chosen so that each patch contains between 50 and 200 cones. Those 52 images were presented in a random order to two ophtalmogist physicians of XV-XX hospital (Kiyoko Nakashima and Michel Paques). They marked the cones on these images thanks to a developped imageJ plugin which enabled us to recover the supposedly true cone positions. The comparisons of the method results with the gold standard were based on (i) the sensibility (Se) index which is the percentage of identifed cones by the algorithm which are well inside the gold standard dataset (Se ∈ [0, 1]) and (ii) the false positive rate (FPR) index which is the ratio of the number of identifed cones by the algorithm which are not inside the gold standard dataset (number of false positives) over the number of cones in the gold standard (F P R > 0 and unbounded). In our procedure, we compared two detection positions (manual and automatic detections) with a tolerance of 3 pixels (it corresponds to the minimal cone diameter 2µm).

For each of the 13 selected images, among their 4 patches, 2 were exactly the same ones. We consider that the less important is the difference (computed in terms of sensibility and false positive rate) between the manual detections obtained with these two same patches, the better is the quality of the image. Table 1 presents the images by increasing quality according to this simple criterion.

Worst image

Best image im01 im02 im03 im04 im05 im06 im07 im08 im09 im10 im11 im12 im13 Not only we proposed to compare our algorithm with this gold standard but we also proposed to compare it with two other existing algorithms. In [START_REF] Li | Automated identification of cone photoreceptors in adaptive optics retinal images[END_REF], the authors proposed a multi-step algorithm. First, they propose to work with an above-threshold version of the image in order to eliminate the dim spots. Afterwards, the regional maxima are detected from a low-pass filtered version of the thresholded image in order to eliminate the noise. Finally, a morphological dilation is used to merge close enough detected spots. The structuring element is chosen according to the minimal cone spacing. In [START_REF] Xue | Photoreceptor counting and montaging of enface retinal images from an adaptive optics fundus camera[END_REF], the same kind of procedure is presented but the regional maxima detection is performed by decreasing greylevel regions in order to favor the brightest detected spots and avoid the threshold step of the method in [START_REF] Li | Automated identification of cone photoreceptors in adaptive optics retinal images[END_REF]. First, we present a graph showing, for each studied method (Xue et al. [START_REF] Xue | Photoreceptor counting and montaging of enface retinal images from an adaptive optics fundus camera[END_REF] in blue, Li et al. [START_REF] Li | Automated identification of cone photoreceptors in adaptive optics retinal images[END_REF], Loquinet al. [START_REF] Loquin | Photoreceptor detection in in-vivo adpative optics images of the retina: towards a simple interactive tool for the physician[END_REF] in black and our method in green) the maximum the difference between Se and FPR for the set of possible parameterization for each method. The number of tested parameters were of 3052 for Xue, 1216 for Li, 7680 (30 standard deviations and 256 tolerances) for the method presented in this article and 256 for Loquin 2011. We took the maximum of Se-FPR since the detection is better for higher value of Se-FPR. This remark leads to a simple interpretation of this graph: globally (except for images 6 and 18) our method is better than the ones proposed in [START_REF] Li | Automated identification of cone photoreceptors in adaptive optics retinal images[END_REF][START_REF] Xue | Photoreceptor counting and montaging of enface retinal images from an adaptive optics fundus camera[END_REF]. We can observe a striking improvement of the detection with our algorithm on the bad quality images. Indeed, for patch images 3, 5, 7, 10, 48 and 49 which are all belonging to two worst images (cf. Table 1), Figure 3 shows that our method particularly outperforms the other ones.

Results

Second, we present the ROC (Receiver Operating Characteristic) cloud and its associated logarithmic fit curve. The ROC curve has FPR in abscissa and Se in ordinate. Such curve is supposed to measure the ability of a decision making support to diagnose. The bigger the area under the curve the more discriminative is the proposed tool. It can be noticed that, according to this ROC analysis, our algorithm is more able to detect cones than the compared ones [START_REF] Li | Automated identification of cone photoreceptors in adaptive optics retinal images[END_REF][START_REF] Xue | Photoreceptor counting and montaging of enface retinal images from an adaptive optics fundus camera[END_REF]. From Figure 4, but in contradiction with Figure 3, we could conclude that the approach without deconvolution (the black curve) is more accurate than the method presented in this article (the black curve). This is due to the fact that the number of tested parameterizations with deconvolution is 30 times bigger than without deconvolution. Therefore, the number of aberrant parameters, thus leading to worst results with deconvolution, is also more important. However this comparison bias due to the number of tested parameters can not be advocated for comparing the Li and Xue approach with the present one. Indeed, the number of tested parameters is higher in our approach than any of the other tested approaches.

Conclusion

In this paper, we proposed a new photoreceptor detection algorithm applied to in-vivo AO images of the retina. This algorithm is an extension of our work in [START_REF] Loquin | Photoreceptor detection in in-vivo adpative optics images of the retina: towards a simple interactive tool for the physician[END_REF] where we add a preliminary deconvolution step. Our method is implemented as an imageJ plugin already in use in a clinical setting at the XV -XX hospital. This use shows its early clinical value before any objective evaluation. However, in this paper, we also present a validation procedure of this algorithm which seems to show that our method is more efficient than the available rare litterature in that field especially for poor quality images. It should be noted that the parameterization of our method is not discussed in this article. A more complete article about our method evoking this sensible point is in preparation. For instance, learning methods of the PSF similar to blind deconvolution are currently tested.

Fig. 1 .

 1 Fig. 1. Usual eye fundus image vs. Adaptive Optics image of the retina on a pathological eye

Fig. 2 .

 2 Fig. 2. Illustration of the deconvolution interest. 1a and 2a: results of our new algorithm on the original images 1 and 2. 1b and 2b: deconvolution results of images 1 and 2

LoquinFig. 3 .

 3 Fig. 3. Maximum of Se-FPR for each image and each compared method

  Loquin et al. with déconvolution Loquin et al. 2011 Li et Roorda 2007 Xue et al. 2007

Fig. 4 .

 4 Fig. 4. ROC cloud and fitted (logrithmic) ROC curve for each compared methods

Table 1 .

 1 Increasing quality image ordering

	10	25	2	17	42	20	32	28	21	23	24	27	51
	3	5	16	1	18	13	52	34	26	37	4	31	8
	48	44	14	9	33	46	29	11	15	30	47	43	6
	49	7	40	35	22	45	36	41	39	12	19	38	50

http://www.imagine-eyes.com/content/view/122/124

This work has been supported by French National Research Agency (ANR) through Tecsan Program (project iPhot n o ANR-09-TECS-009). K. Loquin is currently researcher at LITIS, Rouen and LIRMM, Montpellier, France