Computing solutions of linear Mahler equations
Frédéric Chyzak, Thomas Dreyfus, Philippe Dumas, Marc Mezzarobba
CCOMPUTING SOLUTIONS OF LINEAR MAHLER EQUATIONS
FRÉDÉRIC CHYZAK, THOMAS DREYFUS, PHILIPPE DUMAS,AND MARC MEZZAROBBA
Abstract.
Mahler equations relate evaluations of the same function f atiterated b th powers of the variable. They arise in particular in the study ofautomatic sequences and in the complexity analysis of divide-and-conqueralgorithms. Recently, the problem of solving Mahler equations in closed formhas occurred in connection with number-theoretic questions. A difficulty inthe manipulation of Mahler equations is the exponential blow-up of degreeswhen applying a Mahler operator to a polynomial. In this work, we presentalgorithms for solving linear Mahler equations for series, polynomials, andrational functions, and get polynomial-time complexity under a mild assumption.Incidentally, we develop an algorithm for computing the gcrd of a family oflinear Mahler operators. Introduction
Context.
Our interest in the present work is in computing various classes ofsolutions to linear Mahler equations of the form( eqn ) ‘ r ( x ) y ( x b r ) + · · · + ‘ ( x ) y ( x b ) + ‘ ( x ) y ( x ) = 0 , where ‘ , . . . , ‘ r are given polynomials, r > b ≥ Q between certainof these values come from algebraic relations over ¯ Q ( x ) between the functionsthemselves. This direction was continued by several authors. We refer to Pellarin’sintroduction [19] for a historical and tutorial presentation, and to the referencestherein; see also Nishioka [18] for a textbook.Mahler equations are closely linked with automata theory: the generating seriesof any b -automatic sequence is a Mahler function, that is, a solution of a linearMahler equation; see [10, 11]. Mahler functions also appear in many areas at theinterface of mathematics and computer science, including combinatorics of partitions,enumeration of words, and the analysis of divide-and-conquer algorithms.Very recently, functional relations between Mahler functions have been furtherstudied with a bias to effective tests and procedures [4, 5, 6, 12, 21]. Such studiesmotivate the need for algorithms that solve Mahler equations in various classes offunctions. For instance, testing transcendence of a Mahler series by the criterion Date : April 11, 2018.2010
Mathematics Subject Classification.
Primary: 39A06; Secondary: 33F10, 68W30.This project has received funding from the European Research Council (ERC) under theEuropean Union’s Horizon 2020 research and innovation programme under the Grant AgreementNo 648132. MM was supported in part by ANR grant ANR-14-CE25-0018-01 (FastRelax). a r X i v : . [ c s . S C ] A p r F. CHYZAK, TH. DREYFUS, PH. DUMAS, AND M. MEZZAROBBA of Bell and Coons [6] requires to compute truncations of Mahler series to suitableorders. So does the algorithm by Adamczewski and Faverjon [4, 5] for the explicitcomputation of all linear dependence relations over Q between evaluations of Mahlerfunctions at algebraic numbers. Besides, Mahler functions being either rational ortranscendental—but never algebraic—, solving Mahler equations for their rationalfunctions is another natural approach to testing transcendence, and an alternativeto Bell and Coons’ (see further comments on this in §3.6). Similarly, the hypertran-scendence criterion by Dreyfus, Hardouin, and Roques [12] relies on determining ifcertain Mahler equations possess ramified rational solutions.1.2. Related work.
Mahler equations are a special case of difference equations, inthe sense of functional equations relating iterates of a ring endomorphism σ appliedto the unknown function.Algorithms dealing with difference equations have been widely studied. In par-ticular, the computation of rational solutions of linear difference equations withcoefficients polynomial in the independent variable x is an important basic brickcoming up repeatedly in other algorithms. Algorithms in the cases of the usualshift σ ( x ) = x + 1 and its q -analogue σ ( x ) = qx have been given by Abramov [2, 3]for equations with polynomial coefficients: in both cases, the strategy is to computea denominator bound before changing unknown functions and computing the nu-merator as a polynomial solution of an auxiliary difference equation. Bronstein [8]provides a similar study for difference equations over more general coefficient do-mains; his denominator bound is however stated under a restriction ( unimonomialextensions ) that does not allow for the Mahler operator σ ( x ) = x b .Mahler equations can also be viewed as difference equations in terms of the usualshift σ ( t ) = t + 1 after performing the change of variables t = log b log b x . Thisreduction from Mahler to difference equation, however, does not preserve polynomialcoefficients, which means that neither Abramov’s nor Bronstein’s algorithm can beused in this setting.There has been comparatively little interest in algorithmic aspects specific toMahler equations. To the best of our knowledge, the only systematic study is byDumas in his PhD thesis [13]. In particular, he describes procedures for computingvarious types of solutions of linear Mahler equations [13, Chapter 3]. However,beside a few gaps of effectiveness, that work does not take computational complexityissues into account. To a large extent, the results of the present work can beviewed as refinements of it, with a focus on efficiency and complexity analysis. Morerecently, Bell and Coons [6] give degree bounds that readily translate into algorithmsfor polynomial and rational solutions based on undetermined coefficients. Withregard to series solutions, van der Hoeven [22, §4.5.3] suggests an algorithm thatapplies, under hypotheses, to certain equations of the form ( eqn ) as well as tocertain nonlinear generalizations, and computes the first n terms of a power seriessolution in ˜O( n ) arithmetic operations. At least in the linear case and in analogy tothe case of difference equations, this leaves the open question of an algorithm incomplexity O( n ).1.3. Setting.
Our goal in this article is to present algorithms that compute completesets of polynomial solutions, rational function solutions, truncated power seriessolutions, and truncated Puiseux series solutions of ( eqn ). More precisely, let K bea (computable) subfield of C , and suppose ‘ , . . . , ‘ r ∈ K [ x ]. Denote by K (( x / ∗ )) OMPUTING SOLUTIONS OF LINEAR MAHLER EQUATIONS 3
Kind of solutions Algorithm Complexity K [[ x ]], to order b ν c + 1 Alg. 4 O( rdv + r M( v )) K [[ x ]], to order n , when r = O( d ) Alg. 3 O( rd + nrd ) K [ x ] Alg. 6 ˜O( b − r d + M( d )) K (( x /N )) Alg. 7 ˜O( r N d ( d + n )) K (( x / ∗ )) Alg. 7 ˜O( r b r d ( d + n )) K ( x ), when b = 2 Alg. 9 ˜O( d ) K ( x ), when b ≥ b − r d ) Table 1.
Complexity of the solving algorithms presented in the paper, assuming ‘ = 0. the field S + ∞ n =1 K (( x /n )) of formal Puiseux series with coefficients in K . Let M denote the Mahler operator of radix b , that is the automorphism of K (( x / ∗ )) thatsubstitutes x b for x and reduces to the identity map on K . Writing x again forthe operator of multiplication of a series by x , M and x follow the commutationrule M x = x b M . Equation ( eqn ) then rewrites as Ly = 0 where( opr ) L = ‘ r M r + · · · + ‘ in the algebra generated by M and x . We are interested in the algebraic complexityof computing the kernel of L in each of K [ x ], K ( x ), K [[ x ]], and K (( x / ∗ )).We always assume that ‘ r is nonzero. Except where otherwise noted, we alsoassume ‘ = 0. From a decidability viewpoint, the latter assumption is no loss ofgenerality thanks to the following result [13, Cor. 6, p. 36]. Proposition 1.1.
Given a linear Mahler equation of the form ( eqn ) , one cancompute an equation of the same form, with ‘ = 0 , that has exactly the same formalLaurent series solutions—and therefore, the same polynomial solutions and the samerational-function solutions. Note however that this result does not say anything about the cost of reducingto the case ‘ = 0. We give a complexity bound for this step in §4. As it turnsout, this bound often dominates our complexity estimates for the actual solvingalgorithms. Let us therefore stress that all other complexity results are stated underthe assumption that ‘ is nonzero.For 0 ≤ k ≤ r , we denote by v k ∈ N ∪ { + ∞} and d k ∈ N ∪ {−∞} the valuationand degree of the coefficient ‘ k . Let d ≥ max ≤ k ≤ r d k . Polynomials are implicitlyrepresented in dense form, so that polynomials of degree d in K [ x ] have size d + 1.All complexity estimates are given in terms of arithmetical operations in K , whichwe denote “ops”. The complexity of multiplying two polynomials of degree at most n is denoted by M( n ); we make the standard assumptions that M( n ) = O( n ) andthat n M( n ) /n is nondecreasing.Given two integers or polynomials a and b , we denote their gcd by a ∧ b and theirlcm by a ∨ b ; we use V and W for n ary forms. F. CHYZAK, TH. DREYFUS, PH. DUMAS, AND M. MEZZAROBBA
The following identities are used repeatedly in the text. We gather and repeatthem here for easier reference: ‘ r ( x ) y ( x b r ) + · · · + ‘ ( x ) y ( x b ) + ‘ ( x ) y ( x ) = 0 , ( eqn ) L = ‘ r M r + · · · + ‘ , ( opr ) ν = max k ≥ v − v k b k − , µ = v + ν. ( mu-nu )1.4. General strategy and outline.
The article is organized as follows. In §2,we develop algorithms to compute truncated series solutions of equations of theform ( eqn ). We start with an example that illustrates the structure of the solutionspace and some of the main ideas behind our algorithms (§2.1). Then, we introducea notion of Newton polygons, and use it to prove that the possible valuations (resp.degrees) of the solutions of ( eqn ) in K (( x / ∗ )) (resp. K [ x ]) belong to a finite setthat we make explicit (§2.2). We compute a suitable number of initial coefficientsby solving a linear system (§2.4), then prove that the following ones can be obtainediteratively in linear time, and apply these results to give a procedure that computesa complete set of truncated series solutions (§2.5). Finally, we extend the sameideas to the case of solutions in K [ x ] (§2.6) and in K (( x / ∗ )) (§2.7).The next section, §3, deals with solutions in K ( x ). The general idea is to firstobtain a denominator bound, that is a polynomial q such that Lu = 0 with u ∈ K ( x )implies qu ∈ K [ x ] (§3.1). Based on elementary properties of the action of M onelements of K [ x ] (§3.2), we give several algorithms for computing such bounds(§3.3–§3.4). This reduces the problem to computing a set of polynomial solutionswith certain degree constraints, which can be solved efficiently using the primitivesdeveloped in §2, leading to an algorithm for solving linear Mahler equations in K ( x )(§3.5). We briefly comment on a comparison, in terms of complexity, of Bell andCoons’ transcendence test and the approach by solving the Mahler equation forrational functions (§3.6). The net result is that the new approach is faster.Finally, in §4, we generalize our study to the situation where the coefficient ‘ in ( eqn ) is zero. This makes us develop an unexpected algorithm for computing thegcrd of a family of operators, which we analyze and compare to the more traditionalapproach via Sylvester matrices and subresultants.1.5. Acknowledgment.
The authors are indebted to Alin Bostan for helpful dis-cussions and for pointing us to the work of Grigor’ev [14].2.
Polynomial and series solutions
A worked example.
The aim of this section is to illustrate our solvingstrategy in K [[ x ]] and K (( x / ∗ )) on an example that we treat straightforwardly.In radix b = 3, consider the equation Ly = 0 where(2.1) L = x (1 − x + x )(1 − x − x ) M − (1 − x − x − x − x ) M + x (1 + x )(1 − x − x ) . Assume that y ∈ K (( x / ∗ )) is a solution whose valuation is a rational number v .The valuations of ‘ k M k y , for k = 0 , ,
2, are respectively equal to 6 + v, v, v .If one of these rational numbers was less than the other two, then the valuationof the sum P k =0 ‘ k M k y would be this smaller number, and Ly could not be zero. OMPUTING SOLUTIONS OF LINEAR MAHLER EQUATIONS 5
Consequently, at least two of the three rational numbers 6 + v, v, v have to beequal to their minimum. After solving, we find v ∈ {− / , } .First consider the case v = 3, and write y = P n ≥ y n x n . For m from 10 to 15,extracting the coefficients of x m from both sides of 0 = ‘ y + ‘ M y + ‘ M y , wefind that y , . . . , y satisfy(2.2) 0 = y + y , y + y , − y + y + y , y + y , y + y , − y + y + y . More generally, extracting the coefficient of x m yields the relation(2.3) (cid:0) y m − + y m − − y m − − y m − − y m − − y m − (cid:1) − (cid:0) y m − y m − − y m − − y m − − y m − (cid:1) + (cid:0) y m − − y m − + y m − − y m − − y m − (cid:1) = 0 , where y s is understood to be zero if the rational number s is not a nonnegativeinteger. This equation takes different forms, depending on the residue of m modulo 9:for example, for m = 20 and m = 42, it reduces to, respectively, y + y = 0 , y + y − y − y − y − y − y = 0 . Despite these variations, for any m ≥
10 the index n = m − m ≥
10, we can iterativelyobtain y n from (2.3) in terms of already known coefficients of the series. Conversely,any sequence ( y n ) n ≥ that satisfies (2.3) gives a solution y = P n ≥ y n x n of ( eqn ).As a consequence, the power series solution is entirely determined by the choiceof y and the space of solutions of ( eqn ) in K [[ x ]] has dimension one. A basisconsists of the single series(2.4) x − x + x − x + 2 x − x + 3 x − x + 3 x − x + · · · . The other possible valuation, v = − /
2, is not a natural number. To revert tothe simpler situation of the previous case, we perform the change of variables x = t followed by the change of unknowns y ( t ) = ˜ y ( t ) /t . The equation becomes ˜ L ˜ y = 0with(2.5) ˜ L = (1 − t + t )(1 − t − t ) M − (1 − t − t − t − t ) M + t (1 + t )(1 − t − t ) . To understand this calculation, remember that M was defined on K (( x / ∗ )), so that M ( t ) = M ( x / ) = x / = t .We now expect ˜ L to have solutions ˜ y = P n ≥ ˜ y n t n of valuation 0 and 7 withrespect to t , and the solutions of ˜ L with valuation 0 to correspond to the solutionsof L with valuation − /
2. Extracting the coefficients of x m for m from 0 to 24 from F. CHYZAK, TH. DREYFUS, PH. DUMAS, AND M. MEZZAROBBA both sides of ˜ L ˜ y = 0 and skipping tautologies, we find that ˜ y , . . . , ˜ y satisfy0 = − ˜ y , − ˜ y − ˜ y , y − ˜ y , y − ˜ y , − ˜ y , y + ˜ y , y + ˜ y , y + ˜ y − ˜ y , y + ˜ y , y + ˜ y , y + ˜ y , y + ˜ y , − ˜ y + ˜ y + ˜ y , − ˜ y + ˜ y . Reasoning as above, we derive that, given ˜ y while enforcing ˜ y = 0, there is exactlyone power series solution to ˜ L . More specifically when ˜ y = 1 and ˜ y = 0, we findthe series 1 − t + t − t + t − t + t + · · · . Hence, there is a 2-dimensional solution space in K (( x / ∗ )) for the originalequation ( eqn ), with a basis consisting of the power series (2.4) and the additionalPuiseux series x − / − x / + x / − x / + x / − x / + x / + · · · . Valuations and degrees.
Let us assume that y ∈ K (( x / ∗ )) is a solutionof ( eqn ), whose valuation is a rational number v . The valuation of the term ‘ k M k y is then v k + b k v . Among those expressions, at least two must be minimal to permitthe left-hand side of ( eqn ) to be 0: therefore, there exist distinct indices k , k between 0 and r such that(2.6) v k + b k v = v k + b k v = min ≤ k ≤ r v k + b k v. This necessary condition for Ly = 0 can be interpreted using a Newton polygon analogous to that of algebraic equations [23, Sec. IV.3.2-3]: to each monomial x j M k in L , we associate the point ( b k , j ) in the first quadrant of the Cartesian planeendowed with coordinates U and V (see Fig. 1). We call the collection of these pointsthe Newton diagram of L , and the lower (resp. upper) boundary of its convex hullthe lower (resp. upper ) Newton polygon of L . That two integers k , k satisfy (2.6)exactly means that ( b k , v k ) and ( b k , v k ) belong to an edge E of slope − v of thecorresponding lower Newton polygon.Given an edge E as above, an arithmetic necessary condition holds in addition tothe geometric one just mentioned: the coefficients of the monomials of L associatedto points of E must add up to zero. We call an edge with this property admissible . Example 2.1.
The lower Newton polygon of the operator (2.1) appears in dashedlines in Figure 1. It contains two admissible edges, corresponding to the valuations 3and − / OMPUTING SOLUTIONS OF LINEAR MAHLER EQUATIONS 7 U V µ = 9 P = (1 , v ) Λ : V = − ν U + µ Figure 1.
The Newton diagram of the equation treated in §2.1 for radix b = 3, withcorresponding lower Newton polygon (dashed line) and upper Newton polygon(dotted line). Lemma 2.2.
Let L be defined as in ( opr ) . The valuation v of any formal Puiseuxseries solution of ( eqn ) is the opposite of the slope of an admissible edge of thelower Newton polygon of L . It satisfies − v r b r − ( b − ≤ v = − v k − v k b k − b k ≤ v b − , where ( b k , v k ) and ( b k , v k ) are the endpoints of the implied edge.Proof. The fact that v is the opposite of a slope together with its explicit formfollow from (2.6) and the discussion above. There remains to prove the upper andlower bounds. The leftmost edge of the lower Newton polygon of L provides thelargest valuation and its slope ( v k − v ) / ( b k −
1) for some k ≥ − v / ( b − v r − v k ) / ( b r − b k ) for some k < r , is bounded above by v r / ( b r − b r − ). (cid:3) Proposition 2.3.
The dimension of the space of solutions of the homogeneousequation Ly = 0 in K (( x / ∗ )) is bounded by the order r of L .Proof. The space of solutions admits a basis consisting of Puiseux series withpairwise distinct valuations. The number of possible valuations is bounded by theedge count of the lower Newton polygon of L , which is at most r . (cid:3) F. CHYZAK, TH. DREYFUS, PH. DUMAS, AND M. MEZZAROBBA X Y ∆ = P ∗ Λ ∗ = ( ν, µ ) Figure 2.
The infinite matrix R corresponding to the example treated in §2.1: solidcircles denote nonzero entries, hollow circles denote recombinations to zero. Remark . As we will see, the dimension of the solutions in K (( x / ∗ )) can bestrictly less than r . It is natural to ask how to construct a “full” system of r linearlyindependent formal solutions in some larger extension of K ( x ). We will not pursuethis question here and point to Roques’s work for an answer; see [21, Lemma 20and Thm 35] and [20, Theorem 1]. See also Remark 2.18 below.In analogy with the previous discussion on valuations of solutions, if a Puiseuxseries solution of ( eqn ) involves monomials with maximal exponent δ , then theexpression d k + b k δ must reach its maximum at least twice as k ranges from 0 to r .As we see by the same reasoning as above (or by changing x to 1 /x , which exchangesthe lower and upper Newton polygons), − δ is then one of the slopes of the upperNewton polygon of L . The largest possible value corresponds to the rightmost edge. Lemma 2.5.
The maximum exponent δ of a monomial in a finite Puiseux seriessolution, and in particular the degree of a polynomial solution, is the opposite of theslope of an admissible edge of the upper Newton polygon. It satisfies δ = − d k − d k b k − b k ≤ db r − ( b − , for some k = k . OMPUTING SOLUTIONS OF LINEAR MAHLER EQUATIONS 9
The admissibility of an edge of the upper Newton polygon is defined in analogywith admissibility in the lower Newton polygon.2.3.
The nonhomogeneous case.
One of the proofs of results about Puiseuxseries solutions in §2.7 makes use of extended Newton diagrams that take intoaccount the right-hand side of nonhomogeneous equations.For L as in ( opr ) and a Puiseux series ‘ −∞ of valuation v −∞ ∈ Q ∪ { + ∞} ,consider the nonhomogeneous equation(2.7) ‘ r ( x ) y ( x b r ) + · · · + ‘ ( x ) y ( x b ) + ‘ ( x ) y ( x ) = ‘ −∞ ( x ) . Given a Puiseux series solution y ∈ K (( x / ∗ )) of this equation, with valuation v ∈ Q ,we define the Newton diagram of ( L, ‘ −∞ ) as the Newton diagram of L , augmentedwith all points (0 , α ) for which x α appears with nonzero coefficient in ‘ −∞ . Thenotion of lower Newton polygon extends correspondingly.As in §2.2, these definitions are motivated by analyzing the minimum of thevaluations v k + b k v of the terms of the left-hand side of (2.7): either this minimumis equal to v −∞ , or it is less than v −∞ and must be reached as least twice on theleft-hand side. In both cases, making the convention that b −∞ = 0, there existdistinct indices k , k , now in {−∞ , , , . . . , r } , such that the analogue v k + b k v = v k + b k v = min k ∈{−∞ , , ,...,r } v k + b k v of (2.6) holds. Again, this exactly means that ( b k , v k ) and ( b k , v k ) belong to anedge E of slope − v of the lower Newton polygon, now of ( L, ‘ −∞ ).Depending on v −∞ and ˆ v = min ≤ k ≤ r ( v k + b k v ), the lower Newton polygonof ( L, ‘ −∞ ) can: be equal to that of L , if ‘ −∞ = 0; add an edge to its left,if v −∞ > ˆ v ; prolong its leftmost edge, if v −∞ = ˆ v ; or replace some of its leftmostedges, if v −∞ < ˆ v . We defined the admissibility of an edge E of the lower Newtonpolygon of L in terms of the coefficients of those monomials x v k M k in L associatedto points on E . We extend the definition to edges of the lower Newton polygonof ( L, ‘ −∞ ) by the convention that, if a point has to be considered for k = −∞ ,the corresponding coefficient is the opposite of the coefficient of x v −∞ in ‘ −∞ .Admissibility is again a necessary condition for v to be a possible valuation of asolution of (2.7).2.4. Approximate series solutions.
We now concentrate on the search for powerseries solutions y ( x ) = y + y x + · · · ∈ K [[ x ]] of ( eqn ). Extracting the coefficientof x m in both sides of it yields a linear equation for the coefficients y n . This linearequation can be viewed as a row, denoted R m , of an infinite matrix R = R ( L ).The matrix R consists of overlapping strips with different slopes. We view itsrow and column indices, starting at 0, as continuous variables Y and X with the Y -axis oriented downwards. Each nonzero term ‘ k ( x ) M k then corresponds to matrixentries in the strip b k X + v k ≤ Y ≤ b k X + d k . By definition of v k and d k , theentries lying on the lines Y = b k X + d k and Y = b k X + v k that delimit the strip arenonzero, except maybe at intersection points of such lines (obtained for different k ).Because of our assumption that ‘ is nonzero, the smallest slope is 1, obtained for k = 0.For large Y , the line Y = X + v becomes the topmost one, and each row R m determines a new coefficient y n uniquely, for n = m − v . Thus, the power series solutions are characterized by a finite subsystem of R . In order to state this factmore precisely in Proposition 2.6 below, define( mu-nu ) ν = max k ≥ v − v k b k − , µ = v + ν. In terms of the Newton diagram, ν and µ are, respectively, the opposite of the slopeand the V -intercept of the leftmost edge of the lower Newton polygon. Note that,as we can deduce from the proof of Lemma 2.2, there is no nonzero power seriessolution when ν <
0, which happens if and only if v is a strict minimum of allthe v k over 0 ≤ k ≤ r . Proposition 2.6.
Assume that ν ≥ . A vector ( y , . . . , y b ν c ) is a vector of initialcoefficients of a formal power series solution (2.8) y = y + · · · + y b ν c x b ν c + y b ν c +1 x b ν c +1 + · · · of ( eqn ) if and only if it satisfies the linear system given by the upper left ( b µ c +1) × ( b ν c + 1) submatrix of R . The power series solution (2.8) extending ( y , . . . , y b ν c ) isthen unique.Proof. A series y = y + y x + · · · is a solution if and only if its coefficients satisfythe system ( R m ) m ≥ . Whenever(2.9) v + n < v + b n, . . . , v + n < v r + b r n, the row R v + n of R is the first one with a nonzero entry of index n . It thendetermines y n in terms of y , . . . , y n − . Condition (2.9) is equivalent to n > ν ,hence, for any given ( y n ) ≤ n ≤ ν , there is a unique choice of ( y n ) n>ν satisfying all theequations R m for m > v + ν = µ . As, when (2.9) holds for a given n , the entries ofindex n of R m with m < v + n are zero, the remaining equations ( R m ) ≤ m ≤ µ onlyinvolve the unknowns ( y n ) ≤ n ≤ ν . (cid:3) We note in passing the following corollary, which is the essential argument in theproof of [21, Theorem 22].
Corollary 2.7.
In case the leftmost edge of the lower Newton polygon of L lieson the axis of abscissas and is admissible, Equation ( eqn ) admits a power seriessolution of valuation 0.Proof. We then have ν = µ = 0, so the only condition to check is that the first entryof R is zero. This is equivalent to the edge being admissible. (cid:3) The geometric interpretation of the quantities µ and ν defined by ( mu-nu ) is aspecial case of a general correspondence between the structure of the matrix R andthe Newton diagram of L via the point-line duality of plane projective geometry.The correspondence stems from the fact that a monomial x j M k of L is associatedboth to a point ( b k , j ) in the Newton diagram and, by considering its action onpowers of x , to the entries of R lying on the line Y = b k X + j . More generally,under projective duality, each point ( U, V ) in the plane of the Newton diagramcorresponds to a line Y = U X + V in the plane of the matrix R , while, conversely,the dual of a point ( X, Y ) is the line V = − XU + Y . A line through two points( U , V ) and ( U , V ) corresponds to the intersection of their duals.In particular, the point P = (1 , v ) corresponds to the right boundary ∆ : Y = X + v of the strip of entries of slope 1 in the matrix R (see Figures 1 and 2). Inthe ( U, V )-plane, the line containing the leftmost edge of the lower Newton polygon
OMPUTING SOLUTIONS OF LINEAR MAHLER EQUATIONS 11
Input:
The linear Mahler equation ( eqn ). A transformation φ of theform (2.10). An integer w . A set E = { m , m , . . . } of row indices,with m < m < . . . . Output:
The submatrix R E = ( R m,n ) m ∈ E ≤ n
Matrix R . passes through that point P = ∆ ∗ . This line is Λ : V = − νU + µ and correspondsto the bottommost intersection Λ ∗ = ( ν, µ ) of ∆ with the right boundary of anotherstrip. Below this intersection, the entries of R lying on ∆ are the topmost nonzeroentries of their respective columns, and, at the same time, the rightmost nonzeroentries of their respective rows: as already observed, each row R m then determinesa new y n . Example 2.8.
For the operator L of §2.1, the right boundaries of the stripsassociated to the three terms of L have equations Y = X +6, Y = 3 X , and Y = 9 X +3respectively (dotted lines in Fig. 2). The first two of them meet at Λ ∗ = (3 , Y = X + 6 becomes the rightmost line for Y >
9. For m ≥
10, the row R m reflects the relation (2.3). In particular, the existence of a power series solution isentirely determined by the small linear system that uses the rows R to R andthe unknowns y to y (gray rectangle on Figure 2). Solving the system yields y = y = y = 0 and y arbitrary. We then recover the results of §2.1: the spaceof solutions of ( eqn ) in K [[ x ]] has dimension one and a basis consists of the singleseries (2.4). The V -intercept of the leftmost edge of the lower Newton polygonis µ = 9, and the corresponding slope is − ν = −
3. In this case, it is both the columndimension of the small system and the valuation of the solution. Observe how thebottom right sector depicted in light gray corresponds to the system starting withequations (2.2): as the top left rectangle imposes y = y = y = 0, the dots on theleft of the sector in light gray play no role in the equations.As we will see, in the situation of Proposition 2.6, the coefficients y b ν c +1 to y b ν c + n of y can be computed from y , . . . , y b ν c in O( n ) ops for fixed L . This motivates tocall the truncation to order O( x b ν c +1 ) of a series solution an approximate seriessolution of ( eqn ).2.5. Power series solutions.
Our goal at this point is to describe an algorithmthat computes the formal power series solutions of ( eqn ), truncated to any specifiedorder. We first explain how to compute the entries of the matrix R . It is convenient, for expository reasons, to frame this computation as an individual step that returnsa sparse representation of a submatrix of R corresponding to a subset of the rows.Indeed, in our complexity model dense matrices could not lead to good bounds. Wetherefore define a matrix representation to be row-sparse if iterating over the nonzeroentries of any given row does not require any zero test in K . Then, the algorithmessentially amounts to an explicit expression for the coefficients of recurrencessimilar (2.3), which can as well be computed on the fly.In view of the computation of ramified solutions (§2.7), Algorithms 1 and 2 acceptas input a K -linear transformation φ to be applied to the operator L . In general, φ will take the form(2.10) φ ( x j M k ) = x αb k + βj − γ M k , α, γ ∈ Z , β ∈ N > , β ∧ b = 1 , with α, β, γ chosen such that φ ( L ) has plain (as opposed to Laurent) polynomialcoefficients. The reader only interested in polynomial, rational, and power seriessolutions of L may safely assume φ = id, i.e., α = γ = 0 , β = 1. Lemma 2.9.
Algorithm 1 computes the submatrix R E obtained by taking the first w entries of the rows of R ( φ ( L )) with index m ∈ E in O (cid:0) ( r + d ) | E | (cid:1) ops. Each rowof R E has at most r + 2 d nonzero entries.Proof. Write ˜ L = P rk =0 ˜ ‘ k ( x ) M k = φ ( L ). Recall that the row R m is obtained byextracting the coefficient of x m in the equality ˜ Ly = 0, where y = P n ≥ y n x n . Moreprecisely, R m,n is the coefficient of y n x m in the series˜ Ly = r X k =0 d X j =0 ˜ ‘ k,j x j ∞ X n =0 y n x b k n = ∞ X m =0 ∞ X n =0 (cid:16) X j + b k n = m ˜ ‘ k,j (cid:17) y n x m . The definition of φ translates into ˜ ‘ k,j = 0 when j αb k − γ (mod β ), and otherwise˜ ‘ k,j = ‘ k,j for j = αb k + βj − γ . Therefore, R m,n is equal to the sum of ‘ k,j for( k, j ) satisfying αb k + βj − γ = m − nb k . For fixed m and k , the coefficient ‘ k,j only contributes when βj ≡ m + γ (mod b k ). Its contribution is then to R m,n with n = b − k ( B − βj ) where B = m + γ − αb k , and we are only interested in0 ≤ n < w , i.e., B − b k w < βj ≤ B . Using the assumption that β is coprime with b ,the condition on βj (mod b k ) rewrites as j ≡ j (mod b k ), where j is the integercomputed at step 2b. Therefore, the loop 2c correctly computes the contributionof ˜ ‘ k to the entries of index less than w of the row R m i , and hence the algorithmworks as stated.The only operations in K performed by the algorithm are one addition andpossibly one comparison (to update the sparse structure) at each loop pass overstep 2(c)i. The total number of iterations of the innermost loop for a given i is atmost r X k =0 (cid:24) d k b k (cid:25) ≤ r + bb − d ≤ r + 2 d and bounds the number of nonzero entries in the row of index m i . The complexityin ops follows by summing over i . (cid:3) According to Proposition 2.6, the number of linearly independent power seriessolutions and their valuations are determined by a small upper upper left submatrixof R . As a direct attempt at solving the corresponding linear system would have toohigh a complexity (see Remark 2.11), our approach is to first find a set of candidate OMPUTING SOLUTIONS OF LINEAR MAHLER EQUATIONS 13
Input:
A linear Mahler operator L of order r . A transformation φ of theform (2.10). Integers h, w ∈ N . A set E = { m , . . . , m w − } with m < · · · < m w − < h , such that the submatrix ( R m i ,j ) ≤ i,j A vector ( f , . . . , f σ ) of polynomials of degree less than w .(1) Construct the row-sparse submatrix S E = ( R m i ,j ) ≤ i,j Solutions over prescribed monomial support. solutions, spanning a low-dimensional vector space that contains the approximateseries solutions, and to refine the solving in a second step. Geometrically, theidea to obtain a candidate solution g = g + g x + · · · is to follow the “profile”of R (more precisely, the right boundary of the overlapping strips described in theprevious section), using a single equation R m to try and compute each coefficient g n from g , . . . , g n − . (That is, for each n , we resolutely skip all but one equationssusceptible to determine g n .) By duality, this corresponds to keeping a varying lineof increasing integer slope in contact with the lower Newton polygon, and havingit “pivot” around it. In this process, the only case that potentially leaves a degreeof freedom in the choice of g n is when column n contains a “corner” of the profile,corresponding to an edge of the Newton polygon. As a consequence, it is enoughto construct at most r independent candidates solutions. The second step thenconsists in recombining the candidates in such a way that the equations R m thatwere skipped in the first phase be satisfied.This strategy is made more precise in Algorithm 2, which will then be specializedto power series solutions (and later to other types of solutions) by a suitable choiceof E , h and w . By construction, Algorithm 2 outputs polynomials of degree lessthan w that are solutions of a subsystem of the linear system induced by L . Thesepolynomials need not a priori prolong into actual solutions. Lemma 2.10. Algorithm 2 runs in O( rwd + r w + r M( h )) ops, and returns abasis of the kernel of the linear map induced by φ ( L ) from K [ x ] Assume ˜ ν ≥ . Algorithm 2, called with h = b ˜ µ c + 1 , w = b ˜ ν c + 1 , E = (cid:0) min k (˜ v k + nb k ) (cid:1) ≤ n 1) + 1 = O(˜ v )and w ≤ v / ( b − 1) + 1 = O(˜ v ) in the formula of Lemma 2.10, the total complexityis as announced. (cid:3) Given an approximate series solution, the next terms of the corresponding seriessolutions can be computed efficiently one by one using simple recurrence formulae. OMPUTING SOLUTIONS OF LINEAR MAHLER EQUATIONS 15 Input: A linear Mahler operator L of order r . A transformation φ of theform (2.10). A polynomial ˆ y = y + · · · + y b ˜ ν c x b ˜ ν c such that φ ( L ) ˆ y =O( x b ˜ µ c +1 ), for ˜ ν and ˜ µ defined by (2.11). An integer n . Output: A polynomial y + · · · + y b ˜ ν c + n x b ˜ ν c + n .(1) Use Algorithm 1 with E = {b ˜ µ c + 1 , . . . , b ˜ µ c + n } , h = b ˜ µ c + n + 1,and w = b ˜ ν c + n + 1 to construct a submatrix R E of R .(2) Solve R E ( y , . . . , y b ˜ ν c + n ) T = 0 for y b ˜ ν c +1 , . . . , y b ˜ ν c + n , by forwardsubstitution, starting with the coefficients y , . . . y b ˜ ν c given on input.(3) Return y + · · · + y b ˜ ν c + n x b ˜ ν c + n . Algorithm 3. Prolonging an approximate series solution to any order. Proposition 2.13. Given an approximate series solution ˆ y = y + · · · + y b ˜ ν c x b ˜ ν c of ( eqn ) , Algorithm 3 computes the truncation to the order O( x b ˜ ν c + n ) of the uniquesolution y of ( eqn ) of the form y = ˆ y + O( x b ˜ ν c +1 ) in O(( r + d ) n ) ops.Proof. By Proposition 2.6, the system to be solved at step 2 is compatible. Accordingto the description of R provided above, the submatrix ( R m,n ) m> b ˜ µ c ,n> b ˜ ν c is lowertriangular, with nonzero diagonal coefficients, so that the system can be solved byforward substitution. As explained in §2.4, the output is a truncation of a solutionof φ ( L ). By Lemma 2.9, the cost in ops of step 1 is O(( r + d ) n ), and each row of S contains at most r + 2 d nonzero entries. Therefore, step 2 costs O(( r + d ) n ) ops. (cid:3) Polynomial solutions. Our goal in this subsection is Algorithm 6, whichcomputes a basis of all polynomial solutions. Lemma 2.5 provides us with an upperbound d/ ( b r − b r − ) + 1 = O( d/b r ) for the degree of any polynomial solution. Beforewe take this into account, we provide an algorithm to compute polynomial solutionswith degree bounded by w ≥ 0, which runs in a complexity that is sensitive to w .In the same way as in Proposition 2.12, to obtain candidate polynomial solutions f = f + · · · + f w − x w − , we set f n = 0 for n ≥ w and then compute f n fordecreasing n by “following” the “left profile” of the matrix R (or, dually, the upperNewton polygon). The corresponding specialization of Algorithm 2 is formalized asAlgorithm 5. Proposition 2.14. Assume ν ≥ . Algorithm 2, called with φ = id and (2.12) h = d + ( w − b r + 1 , E = (cid:0) max k ( d k + nb k ) (cid:1) ≤ n ≤ w , returns a basis of the space of polynomial solutions of ( eqn ) of degree less than w .For w = O( d/b r ) , the algorithm runs in ˜O( wd + M( d )) ops.Proof. The proof is similar to that of Proposition 2.12: the extracted submatrixof R is now upper triangular; the zeros on its diagonal correspond to the admissiblenonpositive integer slopes of the upper Newton polygon; the number of such zerosis not more than r . Both preconditions of Algorithm 2 are therefore satisfiedand Lemma 2.10 applies. Additionally, the choice of h in terms of w is such thatdeg( Ly ) < h whenever deg y < w for a polynomial y . So, the basis returned is thatof the kernel of the map induced by L from K [ x ] A linear Mahler operator L of order r . Output: A basis ( f , . . . , f σ ) of approximate series solutions of L .Let µ , ν be as defined by ( mu-nu ). If ν < 0, return (). Otherwise, callAlgorithm 2 with φ = id, h = b µ c + 1 , w = b ν c + 1 , E = (cid:0) min k ( v k + nb k ) (cid:1) ≤ n Approximate series solutions. Input: A linear Mahler operator L of order r . An integer w ∈ N . Output: A basis ( f , . . . , f σ ) of the polynomial solutions of L of degree lessthan w .Let µ , ν be as defined by ( mu-nu ). If ν < 0, return (). Otherwise, callAlgorithm 2 with φ = id, h = max k d k + ( w − b r + 1 , w, E = (cid:0) max k ( d k + nb k ) (cid:1) ≤ n Polynomial solutions of bounded degree. Input: A linear Mahler operator L of order r . Output: A basis ( f , . . . , f σ ) of all polynomial solutions of L .Call Algorithm 5 with w = j max k d k b r − ( b − k + 1 and return the result. Algorithm 6. Basis of polynomial solutions. For the complexity result, the hypothesis on w implies h = O( d ) and r = O(log b d ),so that the conclusion of Lemma 2.10 specializes to ˜O( wd + M( d )) ops. (cid:3) Remark . The loose bound on w , namely w = O( d/b r ), permits in particularto obtain a result when d is not the maximal degree of the ‘ k , but only boundsthem up to a multiplicative constant. In this case, the complexity announced byProposition 2.14 specializes to the same complexity as in Corollary 2.16. This willbe used for the numerators of rational-function solutions in §3.5.By Lemma 2.5, the degree of any polynomial solution is bounded above by δ = d/ ( b r − b r − ) + 1. Specializing Proposition 2.14 to w = b δ c , we obtain abound for the complexity of computing the whole space of polynomial solutions. Corollary 2.16. Assuming ν ≥ , Algorithm 2, called with φ = id , h = 3 d + 1 , w = j db r − ( b − k + 1 , E = (cid:0) max k ( d k + nb k ) (cid:1) ≤ n ≤ w , OMPUTING SOLUTIONS OF LINEAR MAHLER EQUATIONS 17 computes a basis of the polynomial solutions of ( eqn ) in ˜O( d /b r + M( d )) ops.Proof. Observe that the choice for w induces that h , as defined in Algorithm 6,satisfies h ≤ d + 1. The result follows from this fact and w = O( d/b r ). (cid:3) Puiseux series solutions. We now discuss the computation of solutionsof ( eqn ) in K (( x / ∗ )). Even though Proposition 1.1 does not apply, we still assumethat the coefficient ‘ of L is nonzero. There is no loss of generality in doing so: if L = L M w for some w ∈ N , then the Puiseux series solutions of L are exactly the y ( x b − w ) where y ranges over the Puiseux series solutions of L . Additionally, theorder of L is bounded by that of L , so that the complexity estimates depending onit will still hold (and equations of order zero that result from the transformationwhen r = w have no nontrivial solutions).The computation of solutions y ∈ K (( x /N )) with a given ramification index N issimilar to that of power series solutions. In order to compute a full basis of solutionsin K (( x / ∗ )), however, we need a bound on the ramification index necessary toexpress them all. Lemma 2.17, communicated to us by Dreyfus and Roques, andProposition 2.19 below provide constraints on the possible ramification indices. Lemma 2.17. If y ∈ K (( x / ∗ )) is a Puiseux series such that Ly ∈ K (( x /q )) where q is coprime with b , then y ∈ K (( x /q )) for some q coprime with b .Proof. Let q be the smallest positive integer such that y ∈ K (( x /q )). Set g = q ∧ b and q = q /g , so that M y ∈ K (( x b/q )) ⊂ K (( x /q )). The expression y = ‘ − (cid:0) Ly − ( ‘ + · · · + ‘ r M r − ) M y (cid:1) shows that y ∈ K (( x /q )) where q = q q . By minimality of q , we have q = kq for some k ∈ N , which simplifies to q = kg . Since q was assumed to be coprimewith b , this implies g = 1. (cid:3) Remark . Some non-Puiseux formal series solutions of Mahler equations with ‘ = 0 do involve ramifications of order divisible by b : perhaps the simplest example,akin to [9, p. 64] (see also [1]), is y = x /b + x /b + x /b + · · · , which satisfies( M − x b − )( M − y = 0.The following proposition formalizes, as a consequence of Lemma 2.17 and theproperties of Newton polygons discussed in §2.2, that no ramification is neededbeyond those present in the candidate leading terms given by the Newton polygon.Call N the lower Newton polygon of L , and let Q denote the set of denominators q of slopes (written in lowest terms) of admissible edges of N such that q ∧ b = 1. Proposition 2.19. Any Puiseux-series solution y of Ly = 0 belongs to V = P q ∈ Q K (( x /q )) . In particular, the space of solutions of L in K (( x / ∗ )) is containedin K (( x /N )) , where N ≤ b r − denotes the lcm of the elements of Q .Proof. Let y ∈ K (( x / ∗ )) satisfy Ly = 0, and suppose by contradiction that y con-tains a nonzero term of exponent p /q where p ∧ q = 1 and q does not divideany element of Q . Choose p /q minimal with these properties. Write y = y + y where y consists of the terms of y with exponent strictly less than p /q , so that y ∈ V and y has valuation p /q . Then g = Ly belongs to V , so that there exists q ∈ N for which q ∧ b = 1 and g ∈ K (( x /q )). Since Ly = − g , Lemma 2.17 impliesthat y ∈ K (( x /q )) for some q coprime with b . In particular, q is coprime with b . Since p /q is the valuation of a solution of the equation Lz = − g , its opposite s = − p /q is the slope of an admissible edge E of the lower Newton polygon N g of ( L, − g )(see §2.3). On the other hand, because of the definition of Q and the properties q ∧ b = 1 and q Q , the edge E cannot be an edge of N . Therefore, by thedescription in §2.3, g must be nonzero and the edge E must be the leftmost edgeof N g . The valuation of g ∈ V is thus a rational number p /q (not necessarilyin lowest terms) with q ∈ Q , so that in particular q ∧ b = 1. As s is the slopeof E in N g , it is of the form ( q v k − p ) / ( q b k ) for some k ∈ { , . . . , r } . Then, q divides q b k . As it is coprime with b , this implies that q divides q ∈ Q , acontradiction. We have proved that y belongs to V .Next, it is clear that V is contained in K (( x /N )). Finally, letting ( b k i , v i ) denotethe vertices of N (sorted from left to right as i increases), the lcm N satisfies N ≤ Q i ( b k i +1 − k i − < b r , as claimed. (cid:3) Remark . The bound N < b r is tight, as shown by the example of M r − x ,which admits the solution x / ( b r − .In order to obtain an algorithm that computes a basis of the space of Puiseuxseries solutions, there remains to generalize the results of §2.4–2.5 to the case ofsolutions lying in K (( x /N )) where N is given. Motivated by the structure of thespace V described in Proposition 2.19, we do not require here that N be equal tothe lcm of all elements of Q : setting it to the lcm of any subset of these elementsalso makes sense. For the most part, the algorithms searching for power seriessolutions apply mutatis mutandis when the indices m and n are allowed to takenegative and noninteger rational values. Nevertheless, some care is needed in thecomplexity analysis, so we explicitly describe a way to reduce the computation oframified solutions of L to that of power series solutions of an operator ˜ L .Denote x = t β , and consider the change of unknown functions y ( x ) = t α z ( t ), for α ∈ Z and β ∈ N > to be determined. Observe that M t = t b . If y ( x ) is a solutionof Ly = 0, then z ( t ) is annihilated by˜ L = t − γ L t α = t − γ r X k =0 t αb k ‘ k ( t β ) M k = r X k =0 ˜ ‘ k ( t ) M k where γ ∈ Z can be adjusted so that the ˜ ‘ k belong to K [ t ]. We then have ˜ L = φ ( L )where φ is the K -linear map, already introduced in §2.5, that sends x j M k to(2.13) φ ( x j M k ) = t − γ t βj M k t α = t − γ + βj + αb k M k . Viewing monomials x j M k as points in the plane of the Newton diagram, the map φ induces an affine shearing(2.14) [ φ ] : (cid:18) b k j (cid:19) (cid:18) α β (cid:19) (cid:18) b k j (cid:19) + (cid:18) − γ (cid:19) . As in §2.5, denote by ˜ v k and ˜ d k the valuations and degrees of the coefficients of ˜ L ,and by ˜ µ and ˜ ν the quantities defined by ( mu-nu ) with v k replaced by ˜ v k . Lemma 2.21. Fix an edge S of the lower Newton polygon of L , of slope − p/q for (not necessarily coprime) p ∈ Z and q ∈ N . Let c be the V -intercept of the linesupporting S . Set α = p , β = q , and γ = qc in (2.13) . Then:(a) the operator ˜ L = φ ( L ) has polynomial coefficients; OMPUTING SOLUTIONS OF LINEAR MAHLER EQUATIONS 19 UV ••••• •• UV •• •• ••• Figure 3. The transformation in Example 2.22 puts the edge with slope 1 / L (left) onto the U -axis (Newton polygon of ˜ L , right). (b) its Newton diagram is the image of that of ˜ L by [ φ ] , with the edge S beingmapped to a segment of the U -axis;(c) in terms of those of L , the parameters associated to ˜ L satisfy ˜ d k = − qc + pb k + qd k ≥ ˜ v k = − qc + pb k + qv k ≥ , ˜ ν = qν − p ≥ , ˜ µ = q ( µ − c ) ≥ . Proof. Observe that qc is equal to the common value on S of pU + qV . Since theendpoints of S have integer coordinates, this value is an integer, and hence thecoefficients of ˜ L are Laurent polynomials. The transformation [ φ ] of the Newtonplane maps segments of slope s to segments of slope ( α + βs ) / (1 + 0 · s ) = p + qs ,and in particular maps S to a horizontal segment. By the choice of c , that segmentlies on the U -axis. Since q > 0, images by [ φ ] of points above S lie above [ φ ]( S ).As monomials of L correspond to points lying on or above S , their images by φ aremonomials of nonnegative degree. This proves assertion (a). It follows that ˜ L hasa Newton diagram in the sense of our definition, and it is then clear this Newtondiagram is as stated by (b). The expressions of ˜ v k and ˜ d k in (c) are a consequenceof (2.13), using again the positivity of β . Those of ˜ ν and ˜ µ follow. We alreadyobserved that ˜ v k ≥ 0. Finally, − ˜ ν and ˜ µ are, respectively, the slope and V -interceptof the leftmost edge of the lower Newton polygon ˜ N of ˜ L . Since ˜ N has a horizontaledge, ˜ ν and ˜ µ are nonnegative. (cid:3) Example 2.22. Consider again the Mahler operator L in (2.1) treated for b =3 in §2.1. We already observed that the slopes of the Newton polygon of L are − / α = − β = 2, and γ = − 3, to obtain theoperator ˜ L in (2.5). The slopes of the Newton polygon of ˜ L are − Theorem 2.23. Algorithm 7 runs in O( r M( N d ) + rN ( d + ( r + d ) n )) = ˜O( r N d ( d + n )) ops(assuming a softly linear-time polynomial multiplication) and computes the truncationto order O( x n +1 ) of a basis of solutions of ( eqn ) in K (( x /N )) . Input: A linear Mahler operator L as in ( opr ). A ramification index N ∈ N > . A truncation order n ∈ N . Output: A vector (ˆ y , . . . , ˆ y σ ) of truncated Puiseux series.(1) Compute the slope s and V -intercept c of the rightmost admissibleedge of the lower Newton polygon of L with slope in N − Z .(2) Define φ and ˜ L = φ ( L ) according to (2.13), with α = − N s , β = N ,and γ = N c .(3) Call Algorithm 2 on L and φ , with h = b ˜ µ c + 1 , w = b ˜ ν c + 1 , E = (cid:0) min k (˜ v k + nb k ) (cid:1) ≤ n Solving a Mahler equation in K (( x /N )). Proof. The discussion at the beginning of this section shows that z ( x ) ∈ K (( x / ∗ ))is a solution of the operator ˜ L computed at step 2 if and only if y ( x ) = x − s z ( x /N )is a solution of L . By Lemma 2.2 and the choice of s , solutions of L in K (( x /N ))have valuation at least − s , and hence correspond to solutions of ˜ L lying in K [[ x ]].Since the mapping z y is linear and invertible, a basis of solutions of ˜ L in K [[ x ]]provides a basis of solutions of L in K (( x /N )).Let S be the edge of the Newton polygon of L considered at step 1, so thatthe notation of the algorithm agrees with that of Lemma 2.21. Lemma 2.21(c)then provides expressions various parameters associated to ˜ L in terms of s , c , andquantities that can be read off L . Since ˜ ν is nonnegative, Proposition 2.12 appliesand shows that step 3 computes a basis ( f , . . . , f σ ) of the space of approximatesolutions of ˜ L in K [[ x ]] in O( rd ˜ v + r M(˜ v )) ops. Denote by ( z , . . . , z σ ) the basisof power series solutions of ˜ L such that each z i extends f i . Then, according toProposition 2.13, the series ˆ z i computed at step 4 satisfy z i = ˆ z i + O( x N ( s + n )+1 ),and their computation takes O( σ ( r + d )˜ n ) ops. Finally, the truncated Puiseuxseries returned by the algorithm satisfy ˆ y i = x − s ˆ z i ( x /N ), hence are truncations ofelements of a basis of solutions of ˜ L in K (( x /N )).Steps other than 3 and 4 do not perform any operation in K , so that the cost inops of the algorithm is concentrated in those two steps. Let ( b k , v k ) and ( b k , v k )with k < k be the endpoints of S , so that(2.15) qc = pb k + qv k = pb k + qv k . Lemma 2.21(c) gives ˜ v = qv + p − qc . If p ≥ 0, then (2.15) implies qc ≥ p andhence ˜ v ≤ qv ≤ N d . If, now, p < 0, first observe that since b k ≥ b k , we have − pb k ≤ − p ( b k − b k ) = q ( v k − v k ). It follows that − qc = − pb k − qv k ≤ qv k ,whence ˜ v ≤ q ( v + v k ) ≤ N d . In both cases, we have proved that ˜ v = O( N d ). OMPUTING SOLUTIONS OF LINEAR MAHLER EQUATIONS 21 The complexity estimate for step 3 thus rewrites as O( rN d + r M( N d )) ops. As s ≤ d (because all slopes of the Newton polygon are bounded by d in absolute value)and σ ≤ r , that of step 4 becomes O( rN ( r + d )( d + n )) ops. The total running timeis therefore O( r M( N d ) + rN ( d + ( r + d ) n )) ops. (cid:3) Recall that Q denotes the set of denominators q of slopes, written in lowest terms,of admissible edges of N such that q ∧ b = 1. Corollary 2.24. Algorithm 7 with N set to the lcm of elements in Q , returnsthe truncation to order O( x n +1 ) of a basis of solutions of ( eqn ) in K (( x / ∗ )) in ˜O( r b r d ( d + n )) ops, assuming M( k ) = ˜O( k ) .Proof. This follows by combining Proposition 2.19 with Theorem 2.23. (cid:3) Example 2.25. With b = 3, let us consider the order r = 11 Mahler operator L = x − ( x + x ) M + x M − ( x − x ) M + (1 + x − x − x − x ) M − ( x − x − x − x − x ) M − (1 + x + x ) M + ( x + x − x − x ) M − ( x + x − x − x ) M + x M + x M − x M . Its associated parameters are w = 0, v = 568, and a Newton polygon made from fivesegments, all admissible, with slopes − / − 3, 0, 1 / / 5. Except for1458 = 2 · , the denominators are coprime with b = 3 and their lcm is N = 65. Therightmost slope is s = 221 / α = − β = 65, hence γ = − L = t − ( t + t ) M + t M − ( t − t ) M + ( t + t − t − t − t ) M − ( t − t − t − t − t ) M − ( t − t − t ) M + (1 + t − t − t ) M − ( t + t − t − t ) M + t M + t M − M . We want to find a basis of Puiseux solutions for L with a precision O( x n ) where n = 10 . According to Algorithm 7, this leads us to compute a basis of formalseries solutions for ˜ L with a precision O( x ˜ n ) where ˜ n = 65002873. We first applyAlgorithm 4 with ˜ ν = 3888, ˜ µ = 6321121. The computation shows that the space ofsolutions has dimension 2. We extend the solutions to the requested precision byAlgorithm 3 and we obtain a basis of formal series solutions˜ f ( t ) = 1 + t + t + t + t + t + t + t + t + O (cid:0) t (cid:1) , ˜ f ( t ) = t + t + t + t + t + t + t + t + O (cid:0) t (cid:1) . Reversing the change of variable, we find the basis f ( x ) = x − + x + x + x + x + x + x + x + x + O (cid:0) x (cid:1) ,f ( x ) = x + x + x + x + x + x + x + x + O (cid:0) x (cid:1) . These truncated series satisfy Lf = O( x e ), Lf = O( x e ) with e = v + n = 1000568.3. Rational solutions We now turn to the computation of rational function solutions of Mahler equationsof the form ( eqn ). Our algorithm follows a classical pattern: it first computes a denominator bound , that is, a polynomial that the denominator of any (irreducible)rational solution must divide. Then it makes a change of unknown functions andcomputes the possible numerators using the algorithm of §2.6. As is usual withother functional equations, the denominator bound is obtained by analyzing theaction of the operator L on zeros and poles of the functions it is applied to.3.1. Denominator bounds: setting. We will call a rational function p/ ( x ¯ v q ) inlowest terms if it satisfies the following conditions: ¯ v ≥ p, q ∈ K [ x ] are coprimepolynomials; q (0) = 0; and p (0) can be zero only if ¯ v = 0.Consider a rational solution p/ ( x ¯ v q ) of ( eqn ), written in lowest terms. Wealready know from Lemma 2.2 that ¯ v ≤ v r / ( b r − b r − ), so we are left with theproblem of finding a multiple of q .Write T a = W r − i =0 M i a . We will freely use the fact that T ( ab ) | ( T a ) ( T b ) for all a and b . For any j between 0 and r , multiplying the equation ‘ r ( x ) M r y + · · · + ‘ ( x ) M y + ‘ ( x ) y = 0 , by ( M r x ¯ v ) ( M j q ) W i = j M i q and reducing modulo M j q yields(3.1) M j q | x ( b r − b j )¯ v ‘ j ( M j p ) _ i = j M i q. As q is coprime with p and q (0) = 0, Equation (3.1) with j = r implies(3.2) M r q | ‘ r T q. This relation is our starting point for computing a polynomial q ? , depending onlyon ‘ r , such that q | q ? .The algorithm for this task, presented in §3.3, operates with polynomials over K ,but it may be helpful in order to get an intuition to first consider the case K = C .Assume for simplicity that q is squarefree. Equation (3.2) then says that, if α is azero of q , each of its b r th roots is either a b k th root with k < r of some zero of q or a zero of ‘ r . Thus, when α is not a root of unity, its b r th roots are either zerosof ‘ r or roots of lower order of some other zero of q , whose b r th roots then satisfythe same property. (Compare Lemma 3.4 below.) As q has finitely many zeros, thiscannot continue indefinitely, so, in this case, we will eventually find a zero α whose b r th roots are zeros of ‘ r . A difficulty arises when α is a root of unity, but then atmost one of its b th roots can be part of a cycle of the map ζ ζ b (cf. Lemma 3.6), OMPUTING SOLUTIONS OF LINEAR MAHLER EQUATIONS 23 and a closer examination shows that the b − Properties of the Mahler and Gräffe operators. Going back to the gen-eral case, and before making the reasoning sketched above more precise, let us statea few properties of the action of M on polynomials. Besides M , we consider the Gräffe operator defined by G : K [ x ] → K [ x ] , p Res y ( y b − x, p ( y )) . In other words, Gp is the product p ( x /b ) p ( ζx /b ) · · · p ( ζ b − x /b ) for any primitive b th root of unity ζ . While M maps a polynomial p to a polynomial whose complexzeros are the b th roots of the zeros of p , the zeros of Gp are the b th powers of thezeros of p .As a direct consequence of the definitions, M and G act on degrees by:deg M p = b deg p, deg Gp = deg p. Some other elementary properties that will be useful in the sequel are as follows. Lemma 3.1. For any nonzero i ∈ N , the following relations between M and G holdfor all p, q ∈ K [ x ] :(a) G i M i p = p b i ,(b) p | M i G i p ,(c) p | q ⇐⇒ M i p | M i q .Proof. The case i > i = 1 by changing the radix, since M i (resp. G i ) is nothing but the Mahler (resp. Gräffe) operator of radix b i ; so we set i = 1.The assertions (a) and (b) are direct consequences of the definition of G as a resultant.The direct implication in (c) is clear. For the converse, write the Euclidean division q = up + v . If M q = sM p for some s ∈ K [ x ], then ( M u ) ( M p ) + ( M v ) = sM p ,whence M v = 0 since deg M v < deg M p . (cid:3) Lemma 3.2. If p ∈ K [ x ] is monic irreducible and i ∈ N , then G i p = q e for somemonic irreducible q ∈ K [ x ] and e ∈ N . Furthermore, G i p = p if and only if p divides M i p . If this holds for i > , G j p is monic irreducible for any j ∈ N .Proof. To prove the first point, consider the factorization G i p = cq e · · · q e s s of G i p for monic irreducible and pairwise coprime q j and a nonzero c ∈ K . Because ofLemma 3.1(c), the polynomials M i q e , . . . , M i q e s s are pairwise coprime. We have M i G i p = c ( M i q e ) · · · ( M i q e s s ) , and, by Lemma 3.1(b), p | M i q e j j for some j . It follows that G i p | G i M i q e j j = q e j b i j by Lemma 3.1(a), proving the first point.Now if p | M i p , then G i p | p b i , and necessarily there is e ∈ N such that G i p = p e .In fact, e = 1 and G i p = p as G i p and p have the same degree and p is irreducible.Conversely, if G i p = p , then p divides M i p by Lemma 3.1(b).Assume G i p = p for some i > 0. Let j ∈ N and m ∈ N such that mi ≥ j . Then p = G mi p = G mi − j ( G j p ) is monic irreducible, so that G j p is monic irreducibletoo. (cid:3) Lemma 3.3. Let f ∈ K [ x ] be a nonconstant polynomial with f (0) = 0 . If f and itsderivative f are coprime, so are M f and ( M f ) . Proof. Assume f ∧ f = 1. Applying M to a Bézout relation shows that M f ∧ M ( f ) = 1. Now, ( M f ) = bx b − M ( f ), so a common factor s of M f and ( M f ) must divide x . As x cannot divide M f because x (cid:45) f , the only possibility is that s be a constant. (cid:3) The following lemma generalizes the fact that the iterated b th roots of a complexnumber α = 0 are all distinct, except in some cases where α is a root of unity. Lemma 3.4. Let p ∈ K [ x ] be monic and irreducible. For general K , M i p and M j p are coprime for all i > j ≥ if none of the G i p for i ≥ is equal to p . When K = Q , the same conclusion holds if Gp is not equal to p .Proof. We proceed by contraposition, assuming the negation of the common con-clusion: for monic irreducible p , assume M i p ∧ M j p = 1 for some i > j ≥ k = i − j ≥ 1. Lemma 3.1(c) implies that M k p and p are not coprime. Then p divides M k p and Lemma 3.2 implies that G k p = p . This proves the result for gen-eral K . For K = Q , a further consequence is that the map α α b k is a permutationof the roots of p in ¯ Q . Hence, all roots of p satisfy α B = α for some power B = b e of b , with e > 0. This means that p divides x B − x . If p = x , Gp = p ; otherwise, p is a cyclotomic polynomial Φ a with a | b e − 1, so a ∧ b = 1. Applying the formulain [13, Prop. 4 p. 14] yields M Φ a = Q b | b Φ ab , so that p divides M p . Lemma 3.2now implies Gp = p again, completing the proof. (cid:3) Remark . Over a general subfield K ⊂ C , the cyclotomic polynomial Φ a factors asΦ a = Ψ · · · Ψ s and G acts as a cyclic permutation of the Ψ i . See also [13, Chap. 1]for a detailed description of the case a ∧ b = 1.Lemma 3.4 states a result for polynomials p that are not part of a cycle of themap G . As a matter of fact, a related graph whose structure plays a crucial role inwhat follows is that of the map √ G that maps a monic irreducible p to the uniquemonic irreducible q such that Gp is some power of q : we call this map the radical of G , as it ignores the exponent generally introduced by G . An immediate degreeargument shows that the cycles of G are exactly the cycles of √ G , and consist ofmonic irreducible polynomials only.To find a kind of generalization of Lemma 3.4 that applies to polynomials oncycles of √ G , we can always reduce to its hypothesis G i p = p for nonzero i , by“stepping back one step” in the graph of √ G , thus leaving the cycle. Lemma 3.6. Let f ∈ K [ x ] be a nonconstant polynomial with f (0) = 0 . There existsa monic irreducible factor q ∈ K [ x ] of M f such that G k q = q for all nonzero k ∈ N .Proof. Choose a monic irreducible factor p of f and write M p = q · · · q s for monicirreducible q i . By contradiction, assume that for each i , there is some nonzero k i forwhich G k i q i = q i . It follows that for k = k · · · k s and all i , G k q i = q i . Lemma 3.1(a)implies p b = ( Gq ) · · · ( Gq s ), and because of Lemma 3.2, for all i , Gq i is irreducible.Hence, there exist nonzero e i ∈ N such that Gq i = p e i , with b = e + · · · + e s .Therefore, for each i , q i = G k − p e i , so that, as q i is irreducible, e i = 1, and thusall q i are equal to some same monic irreducible ˜ q . It follows that M p = ˜ q b . As p isirreducible, Lemma 3.3 applies to show that M p ∧ ( M p ) = 1, which is impossible.The result follows by setting q = q i for a suitable i . (cid:3) OMPUTING SOLUTIONS OF LINEAR MAHLER EQUATIONS 25Φ a Φ a Φ a Φ a Φ a ... Φ a ... Φ a Φ a ... Φ a Φ a Φ a ... Φ a Φ a Φ a Φ a ... Φ a ... Φ a Φ a ... Φ a Φ a Φ a ... Φ a Φ a Φ a Φ a ... x ... x − x − x − x − x − x + 2 x + 2... x − x + 4 x − x + 4... x + 2 x + 4 x + 2 x + 4... . . . . . . . . .. . . . . . . . . Figure 4. Graph of the radical √ G of the Gräffe operator for b = 6 in Q [ x ]. Here, a is apositive integer, coprime to b . In general, the graph of √ G consists of a looprooted at x (top left), bi-infinite trees (bottom), and cycles between cyclotomicpolynomials with infinite trees rooted at them (top right). Example 3.7. To suggest the graph structures induced by the Mahler and Gräffeoperators, we depict on Figure 4 the graph of the radical √ G . Applying M tosome vertex p in the graph results in the product of all antecedents under themap. For example, M ( x − ) = ( x − x + 2)( x − x + 4)( x + 2 x + 4), and M Φ a = Φ a Φ a Φ a Φ a . In the second example, Φ a appears to the right as aconsequence of it being mapped to itself by G .The depicted case, b = 6, is typical for Q [ x ]. In particular, all cycles have length 1as a consequence of the second part of Lemma 3.4.3.3. Denominator bounds: algorithm. Armed with the previous lemmas, wecan now prove the key result that leads to our main denominator bound. Still,to avoid repetitions in the proof of Proposition 3.10 below, we first state twointermediate lemmas.The following lemma can be expressed more intuitively as follows: for any ˜ f thatis not on a cycle of √ G , any g that appears on the tree rooted at ˜ f of antecedentsunder √ G is also not on a cycle. Lemma 3.8. Let ˜ f ∈ K [ x ] be monic irreducible and satisfy G i ˜ f = ˜ f for all i > .Further, let g ∈ K [ x ] be monic irreducible and divide M j ˜ f for some j ≥ . Then G i g = g for all i > .Proof. Suppose G i g = g for some i ≥ 1. By Lemma 3.2, G j g is monic irreducible,and since G j g | G j M j ˜ f = ˜ f b j , it must be ˜ f . Thus, G i ˜ f = G i + j g = G j g = ˜ f , incontradiction with the definition of ˜ f . (cid:3) Lemma 3.9. Let s ≥ r − and m ≥ be integers, and let f ∈ K [ x ] be monicirreducible, q ∈ K [ x ] be nonconstant, and ‘ ∈ K [ x ] be nonzero, and such that x (cid:45) q , M s f m | M r q | ‘ T q , and M s − i f ∧ q = 1 whenever ≤ i < r . Then M s f m divides ‘ .Proof. Let h k | M s f m for a monic irreducible h ∈ K [ x ] and k > 0, so that h k | ‘ T q . We prove by contradiction that h is coprime with T q : suppose thereexists some i satisfying 0 ≤ i < r such that h divides M i q . Then, G i h divides both G i M i q and G i M s f , which, upon applying Lemma 3.1(a), are equal to powers of q and M s − i f , respectively. This contradicts the coprimality of q and M s − i f . Weconclude that h k | ‘ , and the conclusion follows upon considering all h k | M s f m . (cid:3) The following proposition will be used implicitly as a termination test in Al-gorithm 8: as long as there exists a nonpolynomial rational solution p/q , thenonconstant polynomial u proved to exist contains (potential) factors of q and canbe used to change unknowns in a way that lessens the degree of ‘ r . An interpretationof the structure of the proof is as follows: • If some factor of q appears out of all cycles of √ G , there exists such a factor u with no other factor of q in the tree rooted at u , and this u satisfies M r u | ‘ . • Otherwise, each factor f of q is on a cycle and leads to some antecedent ˜ f under √ G that is on no cycle, for which f divides G ˜ f . Considering allpossible f and taking multiplicities into account, we construct a polynomial u such that M r − u | ‘ and q | Gu . Proposition 3.10. Let ‘ ∈ K [ x ] be a nonzero polynomial and q ∈ K [ x ] be anonconstant polynomial such that x (cid:45) q and M r q | ‘ T q . Then there exists anonconstant u ∈ K [ x ] such that: • either M r u | ‘ , • or M r − u | ‘ and q | Gu .Proof. We consider two cases, the first one being when there exists a monic irre-ducible f dividing q such that G i f = f for all i > 0. In this case, we first prove thatwe can also assume without loss of generality that M j f ∧ q = 1 for all j > 0. Assumethe contrary: that the gcd is nontrivial for at least one j > 0. By Lemma 3.4,the M j f for j ∈ N are pairwise coprime, and since q has finitely many factors, M j f ∧ q = 1 for at most finitely many j . Set j to the maximal possible value and g to a monic irreducible factor of M j f ∧ q . Lemma 3.8 applied to g and ˜ f = f implies that G i g = g for all i > 0, and g can replace f with the added property onthe M j g . At this point, Lemma 3.9 applies with s = r and m = 1, proving that M r f divides ‘ . The proposition is proved in this case by choosing u = f .In the second case, let q = c Q k f m k k be the irreducible factorization of q , fora nonzero constant c and two-by-two distinct monic irreducible f k , and with, foreach k , some i k > G i k f k = f k . Fix any k . Lemma 3.6 provides a monicirreducible factor ˜ f k ∈ K [ x ] of M f k such that G i ˜ f k = ˜ f k for all i > 0. If M i ˜ f k ∧ q was nontrivial for some i ∈ N , this gcd would contain some monic irreducible factor g ,necessarily equal to some f k , and Lemma 3.8 would contradict the existence of i k .So the polynomials M j ˜ f k are coprime with q for all j ∈ N . Upon setting s = r − m = m k , and g = ˜ f k , M s g m = M r − ˜ f m k k | M r f m k k | M r q , and M s − i g = M r − − i ˜ f k is coprime with q for all i satisfying 0 ≤ i < r , so that Lemma 3.9 proves that M r − ˜ f m k k = M s g m divides ‘ . Additionally, Gg = G ˜ f k | GM f k = f bk , so that Gg is apower of f k , hence f k | Gg = G ˜ f k , and next f m k k | G ˜ f m k k . Gathering the results over OMPUTING SOLUTIONS OF LINEAR MAHLER EQUATIONS 27 Input: A linear Mahler equation of the form ( eqn ). Output: A polynomial q ? ∈ K [ x ].(1) Set ‘ := ‘ r , then repeat for k = 1 , , . . . :(a) write ‘ = P b r − i =0 x i M r f i with f i ∈ K [ x ];(b) set u k := V b r − i =0 f i ;(c) set ‘ := ( ‘/M r u k ) W r − i =0 M i u k until deg u k = 0, at which point set t = k − u := V b r − − i =0 f i where ‘ = P b r − − i =0 x i M r − f i .(3) Return u · · · u t ( G ˜ u ). Algorithm 8. Obtain a denominator bound from ‘ r . all k , the ˜ f k are pairwise coprime because the f k are; it follows that all M r − ˜ f m k k divide ‘ and are pairwise coprime, so that, finally, the product u = Q k ˜ f m k k satisfies M r − u | ‘ and q | Gu . (cid:3) Remark . In the first case of the proof, which builds u satisfying M r u | ‘ , it isof interest to compare the construction with that in the case of usual recurrences [2].The obtained u is extremal, in the sense that no other factor of q can be foundin the tree rooted at it, that is to say by iterating √ G backward from it; this isused to compute u from the leading coefficient ‘ of the Mahler operator. In thecase of usual recurrences, the shift operator S (with respect to the variable n )and its inverse S − play roles similar to M and √ G , respectively. In Abramov’salgorithm for denominator bounds, poles are searched for by considering poles thatare extremal in a class α + Z : in particular, a pole β ∈ α + Z with minimal realpart corresponds to a monic irreducible factor u = n − β such that S r u divides theleading coefficient ‘ of the recurrence operator. Corollary 3.12. When d < b r − , Eq. ( eqn ) has no nonconstant rational solution.Proof. With the notation above, Lemma 2.2 implies ¯ v = 0. If a nonconstant q couldsatisfy Eq. (3.2), Proposition 3.10 would apply, inducing the contradiction b r − ≤ deg ‘ r ≤ d . So q is constant, and Lemma 2.5 applies and proves p is constant. (cid:3) Proposition 3.10 forms the basis of Algorithm 8, which repeatedly searches forfactors of the form M r u to “be removed” from ‘ r (while “adding back” other factorsof strictly smaller degree) and accumulates the corresponding u into the denominatorbound. The update of ‘ at step 1c of each loop iteration can be viewed as a changeof unknown functions of the form y = ˜ y/u k in ( eqn ). The search for factors of theform M r u , respectively M r − ˜ u , uses the following property (for radix b r , resp. b r − ). Lemma 3.13. Let f , . . . , f b − , u ∈ K [ x ] . The polynomial ‘ = M f + x M f + · · · + x b − M f b − is divisible by M u if and only if f , . . . , f b − are all divisible by u .Proof. The “if” part is clear. Conversely, fix i < b , and assume that M u | ‘ . Let ω be a primitive b th root of unity. Then, M u = ( M u )( ω j x ) | ‘ ( ω j x ) for all j , hence x − x − x − x − p , = x − x − p , = x − x − p , = x − x − p = x − x − p = x − x − p , = x + x + 2 x − x + 1 p = x + x + 2 x − x + 1... x + 4 x + 17 x − x + 1... x + 76 x + 5777 x − x + 1 x + 76 x + 5777 x − x + 1...... p , = 512 x − x + · · · x − x + · · · p = 2 x − x + 10 x + · · · ... 4 x + 18 x + · · · ... 64 x + 2952 x + · · · x + 2952 x + · · · ... p , = x + 1 p , = x − x + 1 p , = x − x + 1 p = x − x + 1... ...512 x − x − p , = 2 x − p , = 2 x − p = 2 x − x + 2 x + 14 x + 2 x + 1... 64 x + 8 x + 164 x + 8 x + 1... Figure 5. Portion of the graph of the radical √ G of the Gräffe operator used for theresolution in Example 3.14. M u divides b − X j =0 ω − ij ‘ ( ω j x ) = b x i M f i . As M u ∈ K [ x b ] and i < b , this implies M u | M f i , and u | f i by Lemma 3.1(c). (cid:3) Example 3.14. In this example, we let b = 3 and use Algorithm 8 to analyze thepotential poles in rational-function solutions of an operator L = (cid:0) p ( x ) · · · p ( x ) (cid:1) M + · · · , where the p i are polynomials to be found in Figure 5 and the coefficients of M and M will be disclosed below. In the figure and this example, polynomials oflarge size are truncated to their first few monomials, and in most cases, we writethem in factored form, although polynomials are manipulated in expanded form inthe actual algorithm.Following Algorithm 8, we set ‘ = p · · · p . Step 1 is motivated by the firstcase in Proposition 3.10: it strives to solve (3.2) by finding a factor u of q suchthat M u | ‘ . For each i , the only monic irreducible candidate factor of u that can OMPUTING SOLUTIONS OF LINEAR MAHLER EQUATIONS 29 “cover” p i upon application of M is the polynomial p i, in the figure. However, M p i, consists of all factors on the level of p i with same ancestor p i, . So, forexample, M p , = p and p , can be part of u , whereas M p , is a strict multipleof p so that p , cannot be made part of u . As a matter of fact, for k = 1 in theloop, the algorithm finds u = p , p , p , p , at step 1b, after rewriting ‘ in theform ‘ = − M (cid:0) (2 x − x − x + 1)( x − x − x − x − x − x + · · · ) (cid:1) + · · · + x M (cid:0) x − x − x + 1)( x − x − x − x − x − x + · · · ) (cid:1) at step 1a. Step 1c resets ‘ to a polynomial that factors into p , p , p , p , p p p , p , p . Following the same approach for k = 2, a new phenomenon occurs because of theloops in the graph: the candidate factor p , that would “cover” p , appears in itsown tree on the same level as p , , and thus has to be rejected. It follows that thealgorithm finds u = p , p , at step 1b, after rewriting ‘ in the form ‘ = − M (cid:0) ( x − x − x − x − x − x + · · · ) (cid:1) + · · · + x M (cid:0) ( x − x − x − x − x − x + · · · ) (cid:1) at step 1a. Step 1c resets ‘ to a polynomial that factors into p , p , p , p , p , p , p , p , p . Following the same approach for k = 3 leads to u = 1: no further factor u of q exists and helps solving Eq. (3.2) by ensuring M u | ‘ .This leads to step 2, which is motivated by the second case in Proposition 3.10:Eq. (3.2) now implies M q | q ∧ M q , which is solved by finding ˜ u such that M ˜ u | ‘ .A difference to step 1 is that at step 2, candidates are looked for just 2 − u = p , p , p , p , , after rewriting ‘ in the form ‘ = M (cid:0) (2 x − x − x + 1)( x − x − x − x − x − x + · · · ) (cid:1) − xM (cid:0) (2 x − x − x + 1)( x − x − x − x − x − x + · · · ) (cid:1) + x M (cid:0) (2 x − x − x + 1)( x − x − x − x − x − x + · · · ) (cid:1) . From these factors, only p , is cyclotomic. But as the algorithm does not factorpolynomials, the other factors cannot be discarded.At step 3, the algorithm returns the bound q ? = u u G ˜ u = p , p , p , p , p , , where the “+1” indicate factors that could have been saved if a cyclotomic testhad been available. The operator L was indeed constructed so as to admit the twoexplicit rational solutions2 x (2 x − x − x − 1) and x − x − x + 1)( x − x − x − x − , whose denominators are effectively “covered” by q ∗ .We remark that, during the steps of the algorithm, the degree of ‘ has droppedfrom its initial value 145 down to 84, then to 62. Example 3.15. Let b = 3 and let us consider the Mahler equation L = (2 x − x − x + 3)(2 x − x − x − M − ( x + 1)(2 x − x + 1)( x − x − x − x − x + 3) M + x (2 x − x + x + 1)( x − x + 1)( x − x − x − x − x + 3) . Following Algorithm 8, we expand (2 x − x − x + 3)(2 x − x − x − 1) to get ‘ ,which step 1a rewrites ‘ = M (6 x − x − x + 3) + xM ( − x + 3 x + x − 1) + x M ( − x + 3 x + x − 1) + x M (4 x − x − x + 2) . (That is, f = f = f = f = f = 0.) We get u = 2 x − x − x + 1, whichfactors into (2 x − x − x − ‘ to a polynomial that factors into(2 x − x − x − x − x − x − x + 3)( x − x − ‘ as instep 1a, we now find ‘ = M (3 − x ) + xM ( − − x ) + x M ( − − x ) + x M (5 + 40 x ) + x M (5 − x ) + x M (2 x + 9) + x M ( − − x ) + x M (15 + 8 x ) + x M (23) , so that u = 1. We pass to step 2, which expands ‘ in the form ‘ = M ( − x + 40 x − x − x + 5 x + 3) + xM (8 x − x − x + 15 x + 5 x − 4) + x M (2 x − x + 23 x + 9 x − , and ˜ u = (2 x − x − x − q ? = u G ˜ u = (2 x − x − x − x − x − x − y = p/ ( x ¯ v q ) is solution of Ly = 0, where ¯ v ≥ p, q ∈ K [ x ]satisfy x ∧ q = p ∧ q = p ∧ x ¯ v = 1, then q divides q ? . Using the results of §2.2, we findthat 0 could not be a pole of a solution in K ( x ) and therefore ¯ v = 0. Consequently, q ? is a denominator bound. Proposition 3.16. Algorithm 8 runs in O((deg ‘ r ) M( d ) log d ) ops if b = 2 , resp.in O( b − r (deg ‘ r ) M( d ) log d ) ops if b ≥ , and computes a polynomial q ? of degreeat most deg ‘ r if b = 2 , resp. at most (deg ‘ r ) /b r − if b ≥ , such that any rationalfunction solution y of ( eqn ) can be written in the form y = p/ ( x ¯ v q ? ) for some p ∈ K [ x ] and ¯ v ∈ N .Proof. For each k ≥ ‘ k denote the value of ‘ consideredat step 1a, so that the value assigned at step 1c is ˜ ‘ k +1 . (In particular, ˜ ‘ = ‘ r .)First, observe that, after step 1b in each loop iteration, u k is by Lemma 3.13 apolynomial of maximal degree such that M r u k | ˜ ‘ k . In particular, the next value,˜ ‘ k +1 , computed at step 1c, is a polynomial. Set ρ = b r − b r − b − , which is at least 1.Step 1c decreases the degree of ‘ bydeg ˜ ‘ k − deg ˜ ‘ k +1 ≥ deg M r u k − deg T u k ≥ ρ deg u k ≥ ρ. In particular, the loop terminates after at most ρ − (1 + deg ‘ r ) iterations, andtherefore the whole algorithm terminates as well. Second, after step 2, ˜ u is similarly OMPUTING SOLUTIONS OF LINEAR MAHLER EQUATIONS 31 a polynomial of maximal degree such that M r − ˜ u | ˜ ‘ t +1 . Therefore, b r − deg ˜ u isbounded above by the degree of ˜ ‘ t +1 , so that (cid:18) t X k =1 ρ deg u k (cid:19) + b r − deg ˜ u ≤ (cid:18) t X k =1 deg ˜ ‘ k − deg ˜ ‘ k +1 (cid:19) + deg ˜ ‘ t +1 ≤ deg ‘ r , where t denotes, as in Algorithm 8, the last value of k for which deg u k > q ? = u · · · u t ( G ˜ u ). If b = 2, then ρ = 1 anddeg q ? = (cid:0)P k deg u k (cid:1) + deg ˜ u is bounded by deg ‘ r ; if b ≥ 3, then ρ = b r − (cid:18) b − b − b − (cid:19) + 1 b − ≥ b r − and deg q ? is bounded by b − ( r − deg ‘ r .Assume that p/ ( x ¯ v q ) is a solution written in lowest terms. Set ˜ q = q and, for k between 1 and t , define the polynomials ˜ q k = ˜ q k − / ( u k ∧ ˜ q k − ). Let us prove byan induction on k that, for 1 ≤ k ≤ t + 1: (i) x (cid:45) ˜ q k − ; (ii) M r ˜ q k − | ˜ ‘ k T ˜ q k − ;and (iii) q | u · · · u k − ˜ q k − . Initially when k = 1, we have ˜ q = q and ˜ ‘ = ‘ r ,so the three properties hold by our assumption on a solution and Equation (3.2).Assume now that x (cid:45) ˜ q k − , M r ˜ q k − | ˜ ‘ k T ˜ q k − and q | u · · · u k − ˜ q k − . It followsfrom ˜ q k − = ( u k ∧ ˜ q k − ) ˜ q k | u k ˜ q k that x (cid:45) ˜ q k and T ˜ q k − | ( T u k ) ( T ˜ q k ). Furthermore,(3.3) ( M r ( u k ∧ ˜ q k − )) ( M r ˜ q k ) = M r ˜ q k − | ˜ ‘ k T ˜ q k − | ˜ ‘ k ( T u k ) ( T ˜ q k ) . Write M r u k = a k M r (˜ q k − ∧ u k ) and ˜ ‘ k = b k M r u k , for suitable polynomials a k and b k . Upon division by M r ( u k ∧ ˜ q k − ), Equation (3.3) becomes(3.4) M r ˜ q k | a k b k ( T u k ) ( T ˜ q k ) . By construction, a k and M r ˜ q k are coprime, as they are the cofactors of M r ( u k ∧ ˜ q k − )in, respectively, M r u k and M r ˜ q k − , so Equation (3.4) finally becomes M r ˜ q k | ˜ ‘ k M r u k ( T u k ) ( T ˜ q k ) = ˜ ‘ k +1 T ˜ q k . By the divisibility assumption on q and the definition of ˜ q k , q | u · · · u k − ˜ q k − = u · · · u k − ( u k ∧ ˜ q k − ) ˜ q k | u · · · u k ˜ q k , completing the proof by induction.The loop terminates when ‘ no longer has any nonconstant factor of the form M r u ,with ‘ = ˜ ‘ t +1 . At this point, M r ˜ q t | ˜ ‘ t +1 T ˜ q t and q | u · · · u t ˜ q t . If ˜ q t is constant,then q | u · · · u t | q ? . On the other hand, if ˜ q t is not constant, Proposition 3.10applies, as x (cid:45) ˜ q t , which implies that ˜ ‘ t +1 admits a factor of the form M r − u such that ˜ q t | Gu . By Lemma 3.13, step 2 computes a polynomial ˜ u such that M r − u | M r − ˜ u . It follows by Lemma 3.1(c) that u | ˜ u , next that ˜ q t | G ˜ u , so that q divides q ? , again.Let us turn to the complexity analysis. Applying M to a polynomial requiresno arithmetic operation. Each execution of step 1b amounts to b r − d/b r , for a total cost of O(M( d ) log d ) ops.The same argument applies to step 2. Similarly, the chain of lcms at step 1c requiresO (cid:18) r − X i =0 M( b i deg u k ) log( b i deg u k ) (cid:19) = O(M( d ) log d ) ops , as ( P r − i =0 b i ) deg u k = O( d ). Since there are at most ρ − (1 + deg ‘ r ) iterations ofsteps 1b and 1c, the cost of step 1 is O( ρ − (deg ‘ r ) M( d ) log d ). If b = 2, then ρ = 1 and the cost of step 1 is O((deg ‘ r ) M( d ) log d ). If b ≥ 3, then ρ ≥ b r − andthe cost is O( b − r (deg ‘ r ) M( d ) log d ).The computation of G ˜ u from ˜ u at step 3 can be performed in O(M( bd )) ops [7, 15]and the final product can be computed in O(M( d ) log d ) ops using a product tree. (cid:3) Proposition 3.16 implicitly provides a bound on deg q that essentially (when ¯ v = 0and ˜ u = 1, exactly) matches that of Bell and Coons [6, Proposition 2]. However, atighter bound holds, especially for b = 2. Proposition 3.17. With the notation above, q has degree at most ‘ r /b r .Proof. Let g = M r q ∧ T q . On the one hand, (3.2) implies M r q | ‘ r g , so that b r deg q ≤ deg ‘ r + deg g . On the other hand, M g divides h = M r q ∨ T q bydefinition of T , hence g M g divides gh = ( M r q ) ( T q ), whence( b + 1) deg g ≤ b r +1 − b − q. Comparing the two inequalities leads todeg q ≤ ( b − 1) deg ‘ r b r +2 − b r +1 − b r + 1 ≤ ( b − 1) deg ‘ r b r ( b − b − ≤ ‘ r b r since ( b − / ( b − b − ≤ b ≥ (cid:3) Remark . The previous discussion to find q ? is entirely based on (3.1) in thecase j = r and on expressing the solution y with a minimal denominator x ¯ v q . Notingthat (3.1) actually holds also for j = r and even if p ∧ q = 1, we may apply it with0 ≤ j ≤ r − p/ ( x ¯ v q ? ) to get additionalconstraints involving ‘ , . . . , ‘ r − that can be used to remove some factors from q ? .3.4. An alternative bound. We now describe an alternative method for com-puting denominator bounds. While it yields coarser bounds, our estimate for itscomputational cost is better, so that it may be a superior choice in some cases. Theresults of this subsection are not used in the sequel. Proposition 3.19. If x ¯ v q ∈ K [ x ] is the denominator of a rational solution of ( eqn ) written in lowest terms, then it holds that q | ( G r ‘ r ) ( G r +1 ‘ r ) · · · ( G r + K ‘ r ) , K = b log b (3 deg ‘ r ) c − r. Proof. Suppose f is monic irreducible and m is positive such that f m | q , andconsider the condition(3.5) M r + j f | r _ i =0 M i q. Clearly, (3.5) is satisfied for j = 0, while it requires b r + j deg f ≤ b r +1 − b − q, which in turn implies j ≤ log b deg q . Plugging in the bound from Proposition 3.17,we obtain j ≤ log b (3 deg ‘ r ) − r .Choose j maximal such that (3.5) holds. Then M r + j f cannot divide T q , and byLemma 3.3, M r + j f is squarefree. Let h be a monic irreducible factor of M r + j f OMPUTING SOLUTIONS OF LINEAR MAHLER EQUATIONS 33 not dividing T q . In the rest of the proof, we write sqrfree p for the squarefreepart of any polynomial p . For all k ≥ 0, set h k = sqrfree( G k h ), and denoteby m k the multiplicity of h k as a factor of T q . Thus m is zero by definitionof h . Continuing with k ≥ 0, Lemma 3.1(b) implies G k h | M G k +1 h , so that h k | sqrfree( M G k +1 h ) | M sqrfree( G k +1 h ) = M h k +1 . As h m k k | T q , we deduce that h m k +1 k | M h m k +1 k +1 | M T q | T q ∧ M r q , then, by using (3.2), h m k +1 k | ‘ r T q . Thedefinition of m k then yields h δ k k | ‘ r for δ k = max( m k +1 − m k , k to the interval j < k ≤ r + j . Then, by Lemma 3.1(b),(3.6) h k | G k h | G k M r + j f = ( M r + j − k f ) b k , and as h k is squarefree, h k divides M r + j − k f . Since f m | q and 0 ≤ r + j − k < r , h mk divides T q , implying m k ≥ m .By Lemma 3.1(b) and Equation (3.6), G r + j − k h k is f b r + j , so that f divides theformer. Then, f δ k | G r + j − k h δ k k | G r + j − k ‘ r . Forming the product of these bounds for k ranging from 0 to j , we get f m | Q jk =0 G k ‘ r , as m ≤ m j +1 and m = 0. The result follows by considering allpossible ( f, m ) such that f m | q . (cid:3) Proposition 3.20. One can compute a polynomial q ∗ ∈ K [ x ] of degree at most d (log b d − r + 2) and such that q | q ∗ in O(M( d log d ) log d ) ops.Proof. If deg ‘ r < b r − , return 1. This is a valid bound by Corollary 3.12. Otherwise,return the bound from Proposition 3.19. As with the previous bound, the G k ‘ r upto k = r + K = O(log d ) can be computed for a total of O(M( bd ) log d ) ops [7, 15].The product then takes O(M( d log d ) log d ) ops. (cid:3) Computing numerators. In order to obtain a basis of rational solutions y of ( eqn ), it suffices to obtain a bound x ¯ v q ? on denominators as in §3.3, to constructan auxiliary equation corresponding to the change of unknown functions y =˜ y/ ( x ¯ v q ? ), and to search for its polynomial solutions ˜ y . We first note the followingconsequence of Lemma 2.5, already proved by Bell and Coons [6, Prop. 2]. Proposition 3.21. If p, q ∈ K [ x ] , not necessarily coprime, satisfy L ( p/q ) = 0 , then deg p is at most deg q + b d/ ( b r − b r − ) c . The procedure to obtain rational solutions is summarized in Algorithm 9. Proposition 3.22. Algorithm 9 computes a basis of rational solutions of its inputequation. Assuming d ≥ b r − , it runs in ˜O( d M( d ) + 2 r d + M(2 r d )) ops when b = 2 and ˜O( b − r d M( d )) ops when b ≥ . Assuming further M( n ) = ˜O( n ) , it runs in ˜O(2 r d ) = ˜O( d ) ops when b = 2 and in ˜O( b − r d ) ops when b ≥ .Proof. Define δ as in step 1, so that δ ≤ d . If δ < b r − , the algorithm will stop afterstep 2. In this case, Corollary 3.12 states that there are no nonconstant rationalsolution. Therefore, the vector space of rational solutions is K when L (1) = 0 and { } otherwise.Otherwise, the algorithm continues with d ≥ b r − . Assume that y ∈ K ( x ) is arational solution of Ly = 0, and let p = x ¯ v q ? y for q ? and ¯ v computed as in step 3. ByProposition 3.16 combined with Lemma 2.2, p is a polynomial. By Proposition 3.21combined with Lemma 2.5, it has degree at most deg( x ¯ v q ? ) + ¯ v = deg q ? + 2¯ v . Input: A linear Mahler equation of the form ( eqn ). Output: A basis of its space of rational function solutions.(1) Set δ = max deg ‘ k .(2) If δ < b r − : return the basis (1) if L (1) = 0, and the empty basis ()otherwise.(3) Compute q ? using Algorithm 8. Set ¯ v = b δ/ ( b r − b r − ) c .(4) For 0 ≤ k ≤ r , set e k = b bδ/ ( b − c − b k ¯ v and˜ ‘ k = x e k ‘ k Y ≤ i ≤ r, i = k M i q ? . Set ˜ L = ˜ ‘ r M r + · · · + ˜ ‘ .(5) Call Algorithm 5 on the equation ˜ Lp = 0, with w = deg q ? + 2¯ v + 1,to compute a basis ( p , . . . , p σ ) of its polynomial solutions of degreeless than w .(6) Return ( p k / ( x ¯ v q ? )) ≤ k ≤ σ . Algorithm 9. Rational solutions Plugging y = p/ ( x ¯ v q ? ) into Ly = 0 and multiplying the resulting equation by thepolynomial x b bδ/ ( b − c Q ri =0 M i q ? , we see that p satisfies ˜ Lp = 0, where ˜ L is definedas in step 4. As b k ¯ v ≤ bδ/ ( b − 1) for k ≤ r , the e k are nonnegative and the ˜ ‘ k arepolynomials. Thus Algorithm 5 applies and, by Proposition 2.14, p belongs to thespan of the p k computed at step 5 of Algorithm 9. Conversely, for all k , the fraction p k / ( x ¯ v p ? ) is a solution of Ly = 0.After step 2, we have b r = O( d ), that is, r = ˜O(1). By Proposition 3.16, the costof step 3 is ˜O( d M( d )) ops when b = 2 and ˜O( b − r d M( d )) ops when b ≥ 3. Define(3.7) ˜ d = 2 b − b − d + b r +1 − b − q ? = ( O(2 r d ) , b = 2 , O( d ) , b ≥ , where the asymptotic bounds follow from Proposition 3.16. Each polynomial ˜ ‘ k defined at step 4 then satisfiesdeg ˜ ‘ k ≤ e k + δ + b r +1 − b − q ? ≤ ˜ d, so its computation as a product of r + 1 factors can be done in O( r M( ˜ d )) ops.This makes a total of O( r M( ˜ d )) = ˜O(M( ˜ d )) ops to compute the ‘ k ’s. Observe aswell that 1 ≤ w = O( ˜ d/b r ). According to Proposition 2.14, step 5 thus requires˜O( b − r ˜ d + M( ˜ d )) ops, which dominates the cost of step 4. Taking the bounds (3.7)into account, we get that step 5 is dominated by step 3 when b ≥ 3, so that the totalcost is ˜O( d M( d ) + 2 r d + M(2 r d )) ops when b = 2 and ˜O( b − r d M( d )) ops when b ≥ n ) = ˜O( n ), this simplifies to the announced complexityestimates. (cid:3) Example 3.23. We continue Example 3.15. We have seen that the denominatorbound is q ? = (2 x − x − x − x − x − x − y = q ? y , so that OMPUTING SOLUTIONS OF LINEAR MAHLER EQUATIONS 35 Ly = 0 if and only if ˜ L ˜ y = 0, where ˜ L = ˜ ‘ M + ˜ ‘ M + ˜ ‘ for˜ ‘ = (2 x − x − x − x − x − x − × (4 x + 2 x + 1)(2 x − x − x + 3)( x + x + 2 x − x + 1) , ˜ ‘ = − (8 x − x + 1)( x − x − x − x + 1) × ( x − x − x + 2 x + 1)(2 x − x − x + 3)( x + x + 2 x − x + 1) , ˜ ‘ = x (2 x − x + x + 1)( x − x + 1)( x − x − × (4 x + 2 x + 1)(2 x − x + x + 2 x − x + 1) × ( x − x − x + 2 x + 1)( x + x + 2 x − x + 1)(2 x − x − x + 3) . We have to compute the complete set of polynomial solutions of ˜ L ˜ y = 0. Thedegree of ˜ ‘ , ˜ ‘ , ˜ ‘ are respectively 16, 46, 54. Using Lemma 2.5, we find that thedegree of a nonzero polynomial solution is necessarily 4 or 5. Following Algorithm 6,we equate the coefficients on both sides of ˜ L ˜ y = 0 up to degree 54, and we obtainthat ˜ y = ˜ y + · · · + x ˜ y is solution of ˜ L ˜ y = 0 if and only if the vector (˜ y , . . . , ˜ y ) issolution of a system of h = 163 equations. A basis of solutions turns out to consistof (2 x − x − x − x − 1) and ( x − x − x − x − x − Ly = 0 consists of12 x − x − x − . Remark . When Mahler equations are considered in difference Galois theory [12,21], the interest tends to be in base fields on which M acts as an automorphism,such as K (( x / ∗ )) and K ( x / ∗ ) = S + ∞ n =1 K ( x /n ). By combining the strategy ofAlgorithm 9 with Proposition 2.19 about possible ramifications, we obtain analgorithm that computes a basis of solutions of ( eqn ) in K ( x / ∗ ). AssumingM( n ) = ˜O( n ), it runs in ˜O(2 r d ) ops when b = 2 and in ˜O( b r d ) ops when b ≥ ‘ is zero.3.6. Testing transcendence. As was announced in the introduction, solvingMahler equations relates to testing the transcendence of Mahler functions. Inparticular, when computing the rational solutions of a Mahler equation ( eqn ) showsthat there are no nonzero rational solutions, this is a proof that all solutions to ( eqn )are transcendental. We compare in this section the complexity of the transcendencetest by Bell and Coons [6] with that of a test by our rational solving.To this end, we briefly sketch Bell and Coons’ “universal” transcendence test [6]and do a complexity analysis of their approach, using our notation and the samelevel of sophistication with regard to algorithms for subtasks. Define κ = (cid:22) ( b − db r +1 − b r + 1 (cid:23) , κ = (cid:22) d/ ( b − b r − (cid:23) , κ = κ + κ + 1 , B = d + κ b r +1 − b − . Bell and Coons [6, Proposition 2 and Lemma 1] show that any rational solution p/q to ( eqn ) without pole at 0 satisfies deg q ≤ κ , deg p ≤ κ + κ , and that ifa series y ∈ K [[ x ]] solves ( eqn ), then either y − p/q = O( x B +1 ) or y = p/q as series. Then, given y = y + y x + · · · , Bell and Coons consider the matrix M = ( y i + j ) ≤ i ≤ κ, ≤ j ≤ B , whose i th row represents the truncation up to O( x B +1 )of the non-singular part of y/x i . To any nonzero ˜ q in the left kernel of M , theyassociate the polynomial q = ˜ q κ + · · · + ˜ q x κ and find a polynomial p of degree at most κ such that y − p/q = O( x B +1 ), therefore such that y = p/q . This leadsto the equivalence that M is full rank if and only if y is transcendental. Bell andCoons’ test therefore consists in computing the truncation of y up to O( x B + κ +1 ),in forming the matrix M , and in determining if M is of full rank, κ + 1. Onlyconsidering the linear-algebra task, which will dominate the complexity, Bell andCoons’ approach takes O( Bκ ω − ) ops, by the algorithm of Ibarra, Moran andHui [16]. When b = 2, we get κ = O( d ), B = O(2 r d ), and a complexity O(2 r d ω );for b ≥ 3, we get κ = O( d/b r ), B = O( d ), and a complexity O( d ω /b ( ω − r ). In eithercase, the dependency in d is in O( d ω ), being not as good as the ˜O( d ) that can beobtained by Algorithm 9, as Proposition 3.22 justifies.In situations where ( eqn ) has nonzero rational solutions, a given series solution y ∈ K [[ x ]] can easily be tested to be one of them, in O( r ω d ) + ˜O( rd ) ops, becauseonly b ν c + 1 = O( d ) initial coefficients of solutions identify them (see §2.4). So inall cases our Algorithm 9 induces a transcendence test in better complexity withrespect to d than with the approach of [6].4. The case ‘ = 0 and an algorithm for computing gcrd’s In this section, we drop the assumption ‘ = 0. More precisely, we consider a linearMahler equation of the form ( eqn ), with ‘ = · · · = ‘ w − = 0 and ‘ r ‘ w = 0. Wecall the integer w the M -valuation of ( eqn ) and d = max k = w,...,r deg ‘ k its degree .We define the M -valuation and the degree of the corresponding operator ( opr )similarly. The goal of this section is to compute a linear Mahler equation with M -valuation equal to 0, such that the new equation and ( eqn ) have the same set ofseries solutions in K (( x )).The algorithm proposed here, Algorithm 11, can be seen as an improvement overan algorithm given by Dumas in his thesis [13, §3.2.1]. In particular, Algorithm 10,borrowed from [13], performs the subtask of splitting an operator of positive M -valuation into a system of operators of zero M -valuation while preserving thesolution set in K (( x )). Dumas’s algorithm next makes use of the right Euclideanstructure of the algebra M ( K ) of linear Mahler operators with coefficients in K ( x ),and transforms the system into a single, equivalent equation by computing a gcrd(greatest common right divisor) via Euclidean divisions. The problem of thisapproach is that the degree of the obtained equation explodes in the process. Toavoid this, we change the second step of algorithm in [13] so as to reuse Algorithm 10and cancellations of trailing instead of leading coefficients.The splitting process of Algorithm 10 is explained in terms of section maps S i ,each of which maps a polynomial in x and M to a polynomial in x and M , andwhose collection plays the role of a partial inverse for M : for 0 ≤ i < b , let S i bethe K -linear map that sends x j M k +1 to x ( j − i ) /b M k if ( j − i ) /b is an integer andto 0 otherwise. Lemma 4.1. Let L be a linear Mahler operator L of the form ( opr ) and havedegree d and positive M -valuation. Then, whenever ≤ i < b , the section S i ( L ) hasdegree at most d/b . Additionally, L can be reconstructed from its sections by (4.1) L = b − X i =0 x i M S i ( L ) . Proof. The degree bound and relation (4.1) are shown by immediate calculations. (cid:3) OMPUTING SOLUTIONS OF LINEAR MAHLER EQUATIONS 37 Input: A linear Mahler operator L with coefficients in K [ x ]. Output: A set of linear Mahler operators with coefficients in K [ x ] and M -valuation zero.(1) If L = 0, return ∅ .(2) If L has M -valuation 0, return { L } .(3) Return the union of the results of calling the algorithm recursivelyon each section S i ( L ) for 0 ≤ i < b . Algorithm 10. Split of ( opr ). Lemma 4.2. Let L be a linear Mahler operator of the form ( opr ) , with order r , M -valuation w , and degree d . Then, Algorithm 10 returns a set of nonzero linearMahler operators of order at most r − w , M -valuation 0, and degree at most db − w .Proof. This is shown by a straightforward induction on w . (cid:3) Instead of considering usual Euclidean divisions according to decreasing powers,which would compute a gcrd as in [13], we use in Algorithm 11 linear combinationsthat kill constant terms: given two nonzero Mahler operators L and L withcoefficients in K [ x ], M -valuation zero, and coefficient of M respectively c and c ,we write R ( L , L ) for the operator c L − c L , whose coefficient of M is zero.We call this operator the interreduction of L and L and a step of the algorithmthat replaces an operator L by an interreduction R ( L , L ) a reduction step . Lemma 4.3. Let L be a system of Mahler operators. Replacing an element L of L by its sections S ( L ) , . . . , S b − ( L ) does not change the set of solutions of L in K (( x )) . Nor does replacing L by the interreduction R ( L , L ) where L , L aredistinct elements of L .Proof. The second claim is obvious. Regarding the first one (already in [13, §3.2.1]),the decomposition (4.1) shows that any common solution of the S i ( L ) is a solutionof L . If, conversely, y is an unramified solution of L , then the x i M S i ( L ) y , 0 ≤ i < b ,have disjoint support, hence S i ( L ) y = 0 for all i . (cid:3) Here, the degree of R ( L , L ) may well be the sum of the degrees of L and L ,but having generated a multiple of M makes it possible to apply splitting and keepdegrees under control. This leads to Algorithm 11, whose correctness and complexityare given in the following proposition.It is worth mentioning that, in general, the equation ˜ L ( y ) = 0 returned byAlgorithm 11 does not have the same set of solutions in K (( x / ∗ )) as the equation L ( y ) = 0. As an example, let b = 2 and consider L = M − xM . We have ˜ L = 1,and the solution space in K (( x / ∗ )) of ˜ L ( y ) = 0 is { } . On the other hand, thesolution space in K (( x / ∗ )) of L ( y ) = 0 is the K -vector space spanned by x / . Proposition 4.4. The operator L has the same set of solutions in K (( x )) as the op-erator ˜ L returned by Algorithm 11. This operator has order ˜ r ≤ r − w , M -valuation 0,and degree ˜ d ≤ db − w . Furthermore, Algorithm 11 runs in O( rb r M( d/b w )) ops. Input: A nonzero linear Mahler operator L of the form ( opr ), order r , M -valuation w , and degree d . Output: A linear Mahler operator ˜ L of order ˜ r ≤ r − w , M -valuation 0, anddegree ˜ d ≤ db − w .(1) Let L be the result of applying Algorithm 10 to L .(2) While L has at least two elements:(a) choose L with highest order in L , then L from L \ { L } ;(b) compute the result L of applying Algorithm 10 to the interre-duction R ( L , L );(c) replace L by ( L \ { L } ) ∪ L .(3) Return the element ˜ L of the singleton L . Algorithm 11. Normalization to ‘ = 0. Proof. Because L = 0 and by construction of Algorithm 10, the initial set L isnonempty. Next, by construction of Algorithm 11, at any time of a run, L isnonempty and contains only elements of outputs from Algorithm 10, so that, byLemma 4.2, if Algorithm 11 terminates, its output must be nonzero and of M -valuation zero. Lemma 4.3 implies that the original operator L , the system L atany time of the run, and therefore the final operator ˜ L , all share the same set ofsolutions in K (( x )).Let us prove the bound on the order and the degree of ˜ L . By Lemma 4.2, theset L computed at step 1 consists of Mahler operators with orders bounded by r − w and degrees bounded by db − w . These bounds keep on holding after each run ofthe loop body at step 2: As the operators L and L chosen at step 2a satisfythe property, their combination R ( L , L ) (including the case it is zero) has orderbounded by r − w , degree bounded by 2 db − w , and positive valuation. By Lemma 4.2,the set L computed at step 2b consists of Mahler operators with orders boundedby r − w − db − ( w +1) . As 2 /b ≤ 1, the set L retainsthe property after the update at step 2c. Therefore, if the algorithm terminates, itreturns at step 3 an element of L , therefore with the announced order and degreebounds.We finally prove termination and complexity by a joint argument. To this end, werepresent the process of Algorithm 11 by an oriented tree labeled by operators L nw ,for integers n and words w on the alphabet { , . . . , b − } . These operators L nw will bethe operators considered during the execution of the algorithm. This tree is rootedat the node labeled L (cid:15) = L , and evolves by following the execution of Algorithm 11.Each time a section of an operator L nw is computed by the subtask of Algorithm 10,whether it be at step 1 or at step 2b, the tree is augmented by new edges from L nw toits subsection L nwj = S j ( L nw ). For each choice of L = L nw and L = L n w at step 2a,the tree is augmented by a new edge labeled L n w , from L nw to L m +1 (cid:15) , if m is thelarger upper index in the tree before reduction. Thus, one obtains that at each stageof the execution, the set L is equal to the collection of nonzero leaves of the currenttree. Now, by construction of the tree and by design of the algorithm, a reductionstep results either in a zero operator or in an operator with positive M -valuation OMPUTING SOLUTIONS OF LINEAR MAHLER EQUATIONS 39 that is immediately split to its sections. Therefore, following a path from the rootto a leaf, two reduction edges can only appear if separated by at least one sectionedge. As section edges reduce orders by at least 1, while reduction edges do notincrease orders, the tree has to be finite and the algorithm terminates. The onlyarithmetic operations of the algorithm are the polynomial products involved in thecomputation of the R ( L , L ) at step 2b. It was proved above that any operatorof L has degree bounded by d/b w . Because operators all have order at most r andas the size of the tree bounds the number of reductions, the algorithm has totalcomplexity O( rb r M( d/b w )). (cid:3) Remark . A slightly better complexity can be obtained by a variant of Algo-rithm 11, in which the L at step 2a is not chosen as having maximal order, butaccording to a notion of depth in the tree introduced for the proof of Proposition 4.4.Doing so guarantees a better behavior of degrees, with a geometric decrease withdepth, as opposed to the uniform bound d/b w used in the proof above.Define the depth β of a node L nw in the tree as the number of section edgesfrom the root L (cid:15) to L nw , and change the strategy at step 2a to choose L amongthe elements of L of lowest depth. By another induction, L nw has order not morethan r − β , as in the proof above, but its degree is not more than d/b w if β ≤ w ,and not more than (2 /b ) β − w ( d/b w ) if β > w . A bound on the complexity becomes r X β = w ( r + 1) b β M (cid:18) β − w db β (cid:19) ≤ O ( r M(2 r d/b w )) . This bound is better than the original complexity O( rb r M( d/b w )) when b ≥ b = 2, the new bound is not tight and the variant algorithm has the samecomplexity bound as Algorithm 11. Example 4.6. We apply Algorithm 11 with b = 3 and the operator L = ‘ M + ‘ M + ‘ M + ‘ M with ‘ = x (1 − x + x + x − x + x )(1 − x + x ) ,‘ = − x (cid:0) x − x − x + x + x + x + x − x − x + x + x + x + x − x − x + x + x + x + x − x − x + x + x + x + x − x − x + x + x (cid:1) ,‘ = (cid:0) x − x + x + x + x − x − x + x + x + x − x + x + x + x + x + x + x − x + x + x + x − x − x + x + x + x − x + x + x (cid:1) ,‘ = − (1 + x + x )(1 − x + x )(1 − x + x + x − x + x ) . Starting from L (cid:15) = L , we compute its sections (see Fig. 6, blue edges): first, L = S ( L (cid:15) ), which has M -valuation 0 so that the process of splitting stopsfor it; next, L = S ( L (cid:15) ), which is zero and is dropped; last, L = S ( L (cid:15) ), whichhas M -valuation 1. Splitting continues for the latter and provides L = S ( L ), L = S ( L ), L = S ( L ), all with M -valuation 0. Note that during this splitting,the operators L (cid:15) , L = 0, and L disappear. A reduction is made (see Fig. 6, red L (cid:15) = L (4 , L (3 , L = 0 L (3 , L (2 , L (2 , L (2 , L (cid:15) (3 , L L (2 , L (2 , L (2 , L (cid:15) = 0 L L (cid:15) = 0 L L (cid:15) = 0 L L (cid:15) = 0 L L (cid:15) = 0 L Figure 6. Execution of Algorithm 11 on the operator of Example 4.6. Each nonzerooperator is given with a corresponding pair (order, degree). Operators aregenerated in the following order: L (cid:15) = L , L , L , L , L , L , L , L (cid:15) , L , L , L , L (cid:15) , L (cid:15) , L (cid:15) , L (cid:15) , L (cid:15) . Blue and red arrows respectively represent sectionand reduction steps. Labels on (red) arrows provide the auxiliary operatorsused for reduction. The process starts with L (cid:15) = L and ends with L . Observethe strict decrease of orders along blue edges and large decrease along rededges. Also observe that degrees are divided by at least 3 on blue edges and,for the only nontrivial red edge of this example, how the reduction of L by L induces an increase of the degree from 49 to 58, which is not morethan 49 + 12. edges) where R ( L , L ) = L (cid:15) replaces L . The process continues and, at the end,there only remains L = x (1 + x + x )(1 − x + x )(1 − x + x ) − x (1 + x + x )(1 − x + x )(1 − x + x − x + x )(1 + 2 x + x ) M + x (1 + x + x )(1 − x + x )(1 + x + x )(1 − x + x ) M . It is worth noting that L has a content c = x (1 + x + x )(1 − x + x ), sothat we can write L = c ¯ L where ¯ L is a primitive polynomial (with respectto M ). The computation shows that L is in the left ideal generated by ¯ L inthe algebra M ( Q ). This and exhibiting the M -valuation w = 1 of L providesfactorizations L = L M = L M ¯ L . We can say that M w has been pushed asmuch as possible to the left. Using Algorithm 9, we find that a basis of solutionsof L in K ( x ) is given by 1 and xx − . Since L has order two, this also formsa basis of solutions of L in K (( x )), as a consequence of Proposition 2.3, and byProposition 4.4, a basis of solutions of L in K (( x )).We now proceed to prove that Algorithm 11 indeed computes a gcrd withcontrolled degree. This is proved in Theorem 4.9 below, using the following lemmas. Lemma 4.7. For any operators P , P , and any integer i such that ≤ i < b , S i ( P M P ) = S i ( P M ) P . OMPUTING SOLUTIONS OF LINEAR MAHLER EQUATIONS 41 Proof. By linearity, it is sufficient to consider P = x j M k and P = x j M k . Then, P M P = x j + b k j M k + k +1 . Either b divides j − i and S i ( P M P ) = x ( j − i ) /b + b k j M k + k = x ( j − i ) /b M k x j M k = S i ( P M ) P , or b does not divide j − i and both extreme terms are zero, thus equal again. (cid:3) Lemma 4.8. For any operators P , P , and P , all of M -valuation 0, let c be thecoefficient of M in P . Then, R ( P P, P P ) = cR ( P , P ) P .Proof. The property holds, as obviously the coefficient of M in a product is theproduct of the coefficients of M in the factors. (cid:3) Theorem 4.9. Steps 2 and 3 of Algorithm 11 compute a gcrd of the elements ofthe split L of L obtained at step 1. The degree of this particular gcrd is bounded bythe maximal degree of the elements of L .Proof. Let I denote the left ideal M ( K ) L generated by L at any time in the runof the algorithm. Call G the monic gcrd of the elements of the set L as obtainedfrom L at the end of step 1. By (4.1), G is a right factor of L . By the definitionof R ( · , · ) and because of (4.1) again, the ideal I can only increase during the run ofthe algorithm, so that during step 2, M ( K ) L ⊂ M ( K ) G ⊂ I .We show by induction that G is a right factor of all elements of L at any time instep 2, in other words, that I ⊂ M ( K ) G . This is true by the definition of G whenentering the loop. The set L contains only elements with M -valuation 0, and itcannot be empty when entering the loop, so G has M -valuation 0 as well. At anystep 2b, divisibility on the right by G is preserved for R ( L , L ), by Lemma 4.8. As R ( L , L ) has positive M -valuation, one can choose P = G and find P so as towrite R ( L , L ) = P M P . By Lemma 4.7, it follows that divisibility on the rightby G is also preserved for each element of L , then for each element of the nextvalue of L .As a consequence, during step 2, I constantly equals M ( K ) G . In particular, thefinal operator ˜ L is proportional to G .The degree bound was proved as part of Proposition 4.4. (cid:3) Note that the origin of the initial L as a split of L , at step 1 of Algorithm 11plays no role in the proof of Theorem 4.9. Thus, Algorithm 11 implicitly contains analgorithm for computing the gcrd of any family L of operators of M -valuation zero. Remark . We developed Algorithm 11 without targeting a gcrd and realizedTheorem 4.9 only a posteriori. As Algorithm 11 indeed works by computing agcrd as the original algorithm in [13], it is now instructive to compare the resultof Proposition 4.4 with bounds on the size of gcrds of Mahler operators given byexisting methods. Such a bound can be computed using a variant of the subresultantargument given by Grigor’ev [14, §5] in the differential case.Let L , . . . , L n be operators of respective order r ≥ r ≥ · · · ≥ r n ≥ d , . . . , d n ≤ δ . Let G = U L + · · · + U n L n be their greatest common right divisor.We can assume that the order of each term U i L i is less than t = r + r n . Indeed, forall i, j the linear equation V i,j L i = V j,i L j with V i,j , resp. V j,i , constrained to havedegree at most r j , resp. at most r i , has nontrivial solutions. Via Euclidean divisions U i = Q i V i,n + R i , we obtain G = P i ( Q i V i,n + R i ) L i = P i Q i V n,i L n + P i R i L i = P i ˜ U i L i where the ˜ U i for i ≤ n − r n . The n − U i L i as well as G itself have order less than r + r n , hence the same must be true of ˜ U n L n . Consider a Sylvester-like matrix S ∈ K [ x ] s × t with rows R ( L ) , R ( M L ) , . . . , R ( M t − r − L ) , . . . , R ( L n ) , R ( M L n ) , . . . , R ( M t − r n − L n ) , where, for any operator L = P k ‘ k M k , we denote R ( L ) = ( ‘ t − , . . . , ‘ ). Call C , C , . . . , C t − the columns of S , listed from right to left (so that C j containsthe coefficients of M j in M k L i ), and C j, , C j, , . . . , C j,s − the entries of C j . Let m denote the order of G , and choose J ⊆ { m + 1 , . . . , t − } of cardinality | J | =rk S − C j with j ∈ J form a basis of the spanof C m +1 , . . . , C t − , while the C j for j ∈ { m } ∪ J form a basis of the full columnspace of S . To see that such a J exists, consider a row echelon form of S : since R ( G ) belongs to the left image of S and G has minimal order among the nonzeroelements of the ideal P i M ( K ) L i , the rightmost pivot lies on column m . Further,let I ⊆ { , . . . , s − } be such that the submatrix ( C j,i ), i ∈ I , j ∈ J ∪ { m } of S isnonsingular. Call D m the corresponding minor, and more generally define D k asthe determinant of the submatrix ( C j,i ), i ∈ I , j ∈ J , extended on the right by acopy of C k . Expanding D k along the last column yields D k = P s − i =0 u i C k,i , wherethe u i do not depend on k . For each k > m , the determinant D k is zero, as C k isin the span of the C j for j ∈ J . It follows that the vector (0 , . . . , , D m , . . . , D )belongs to the left image of S . Thus, there is a gcrd of L , . . . , L n with polynomialcoefficients whose coefficients are minors of S .The entries of S have degree bounded by δ = max ni =1 ( b t − r i − d i ). Therefore,the degree of G is as most tδ ≤ r b r − δ ≤ r b r δ . Using fast polynomial linearalgebra, it is plausible that one could actually compute G based on this approachwith a complexity of the type ˜O( δ t ω ) = ˜O( b r δ ). Now, the gcrd in the algorithmof [13] is that of a family of iterated sections of the input operator L . In terms ofthe order r and degree d of L , this family can involve elements simultaneously oforder r − d/ 2. Thus, Grigor’ev’s approach (at least in a straightforwardway) would lead to a complexity bound similar to that of Proposition 4.4, but anexponentially worse bound on the degree of the output for large r .This result leaves open the question of devising algorithms for computing solutionsof linear Mahler equations that run in polynomial time in r and d , for all possiblecombinations of these parameters, even when the trailing coefficient ‘ of theequation is zero. In particular, it would be interesting to see if the bounds on thesize of an operator equivalent to L implied by Algorithm 10 would be enough toextend the algorithms of §2–3 to the case where ‘ is zero, without going throughthe explicit computation of such an operator.We end the section by providing an extension of Algorithm 11, which computesa gcrd for a family of operators of arbitrary M -valuations. Theorem 4.11. Algorithm 12 computes a gcrd of the input operators L , . . . , L s .Proof. Observe that the minimal M -valuation of operators in a family is the minimal M -valuation of elements of the left ideal generated by the family, in particular,the M -valuation of any gcrd of the family. This justifies the general design of thealgorithm, with the factorization of M w on the right at step 1.By construction, the L i ’s thus obtained have orders at most r − w and degreesat most d , and at least one, say L , has M -valuation zero. Let G denote themonic gcrd of the L i , which, as L , has M -valuation zero. By Lemma 4.7, G is aright-hand factor of all elements of the set L computed at step 2. By a proof similar OMPUTING SOLUTIONS OF LINEAR MAHLER EQUATIONS 43 Input: A finite family { L i } si =1 of linear Mahler operators with polynomialcoefficients, orders at most r , minimal M -valuation w , and degreesat most d . Output: A linear Mahler operator ˜ L of order ˜ r ≤ r , M -valuation w , anddegree ˜ d ≤ d .(1) Write each L i in the form L i M w , for a polynomial L i in x and M .(2) Let L be the union of the results of applying Algorithm 10 tothe L i ’s.(3) While L has at least two elements:(a) choose L with highest order in L , then L from L \ { L } ;(b) compute the result L of applying Algorithm 10 to the interre-duction R ( L , L );(c) replace L by ( L \ { L } ) ∪ L .(4) Write ˜ L for the single element of the singleton L and return ˜ LM w . Algorithm 12. Computation of a gcrd of an arbitrary family. to the one for Theorem 4.9, it remains so for all subsequent values of L , so for the ˜ L of step 4 as well.As ˜ L is also obviously a right-hand factor of all previously computed operators,including the L i ’s, ˜ L is a gcrd of the latter. This concludes the proof. (cid:3) References 1. Shreeram S. Abhyankar, Two notes on formal power series , Proc. Amer. Math. Soc. (1956),no. 5, 903–905.2. S. A. Abramov, Rational solutions of linear differential and difference equations with polynomialcoefficients , Zh. Vychisl. Mat. i Mat. Fiz. (1989), no. 11, 1611–1620, 1757.3. , Rational solutions of linear difference and q -difference equations with polynomialcoefficients , Programmirovanie (1995), no. 6, 3–11.4. Boris Adamczewski and Colin Faverjon, Méthode de Mahler, transcendance et relationslinéaires : aspects effectifs , J. Théor. Nombres Bordeaux (2016), 17 pages. To appear. https://arxiv.org/pdf/1610.09136.pdf .5. , Méthode de Mahler : relations linéaires, transcendance et applications aux nombresautomatiques , Proc. Lond. Math. Soc. (3) (2017), no. 1, 55–90. MR 36699336. Jason P. Bell and Michael Coons, Transcendence tests for Mahler functions , Proc. Amer. Math.Soc. (2017), no. 3, 1061–1070.7. Alin Bostan, Philippe Flajolet, Bruno Salvy, and Éric Schost, Fast computation of specialresultants , J. Symbolic Comput. (2006), no. 1, 1–29.8. Manuel Bronstein, On solutions of linear ordinary difference equations in their coefficientfield , J. Symbolic Comput. (2000), no. 6, 841–877.9. Claude Chevalley, Introduction to the theory of algebraic functions of one variable , Mathemat-ical surveys and monographs, American Mathematical Society, Providence, R.I, 1951.10. Gilles Christol, Ensembles presque periodiques k -reconnaissables , Theoret. Comput. Sci. (1979), no. 1, 141–145.11. Gilles Christol, Teturo Kamae, Michel Mendès France, and Gérard Rauzy, Suites algébriques,automates et substitutions , Bull. Soc. Math. France (1980), no. 4, 401–419.12. Thomas Dreyfus, Charlotte Hardouin, and Julien Roques, Hypertranscendance of solutionsof Mahler equations , J. Eur. Math. Soc. (2015), 26 pages. To appear. http://arxiv.org/abs/1507.03361 . 13. Philippe Dumas, Récurrences mahlériennes, suites automatiques, études asymptotiques , Thèsede doctorat, Université Bordeaux I, 1993, 241 pages. https://tel.archives-ouvertes.fr/tel-00614660 .14. D. Yu. Grigor’ev, Complexity of factoring and calculating the gcd of linear ordinary differentialoperators , J. Symbolic Comput. (1990), no. 1, 7–37.15. Peter Henrici, Applied and computational complex analysis , vol. III, Wiley Interscience, 1986.16. Oscar H. Ibarra, Shlomo Moran, and Roger Hui, A generalization of the fast LUP matrixdecomposition algorithm and applications , J. Algorithms (1982), no. 1, 45–56.17. Kurt Mahler, Arithmetische Eigenschaften der Lösungen einer Klasse von Funktionalgleichun-gen , Math. Ann. (1929), no. 1, 532.18. Kumiko Nishioka, Mahler functions and transcendence , Lecture Notes in Mathematics, vol.1631, Springer Verlag, Berlin, 1996.19. Federico Pellarin, An introduction to Mahler’s method for transcendence and algebraic inde-pendence , t -motives: Hodge structures, transcendence and other motivic aspects (G. Boeckle,D. Goss, U. Hartl, and M. Papanikolas, eds.), European Mathematical Society, 2016, 47 pages.To appear. Proceedings of BIRS, Banff, Canada, 2009. http://arxiv.org/abs/1005.1216 .20. Julien Roques, On the local structure of Mahler modules , , 2016.21. , On the algebraic relations between Mahler functions , Trans. Amer. Math. Soc. (2018), no. 1, 321–355.22. Joris van der Hoeven, Relax, but don’t be too lazy , J. Symbolic Comput. (2002), 479–542.23. Robert John Walker, Algebraic curves , Springer Verlag, New York, Paris, 1978, reprint fromthe 1950 edition, published by Princeton University Press. Frédéric Chyzak, INRIA, Université Paris-Saclay (France) E-mail address : [email protected] Thomas Dreyfus, CNRS (France), Institut de Recherche Mathématique Avancée, UMR7501, Université de Strasbourg, 7 rue René Descartes 67084 Strasbourg E-mail address : [email protected] Philippe Dumas, INRIA, Université Paris-Saclay (France) E-mail address : [email protected] Marc Mezzarobba, Sorbonne Université, CNRS, Laboratoire d’Informatique de Paris 6,LIP6, F-75005 Paris, France E-mail address ::