Table Of Content

Finding Hyperexponential Solutions of Linear ODEs by Numerical Evaluation Fredrik Johansson∗ Manuel Kauers∗ Marc Mezzarobba RISC RISC Inria,Univ. Lyon,AriC,LIP† JohannesKeplerUniversity JohannesKeplerUniversity ENSdeLyon,46alléed’Italie 4040Linz,Austria 4040Linz,Austria 69364LyonCedex07,France [email protected] [email protected] [email protected] 3 1 0 ABSTRACT calledhyperexponentialifthequotientD(y)/ycanbeidenti- 2 We present a new algorithm for computing hyperexponen- fiedwith arational function. Typicalexamples arerational functions(e.g.(5x+3)/(3x+5)),radicals(e.g.√x+1),ex- n tial solutions of ordinary linear differential equations with ponentials (e.g. exp(3x2 4) or exp(1/x)), or combinations a polynomial coefficients. The algorithm relies on interpret- of these (e.g. √x+1exp−(x9/(x 1))). Equivalently, y is J ingformal seriessolutions atthesingular pointsasanalytic − 1 functionsandevaluatingthemnumericallyatsomecommon called hyperexponentialif thereissome first orderoperator q D+q with q ,q polynomials which maps y to zero. If 1 ordinary point. The numerical data is used to determine a 1 0 0 1 we regard differential operators as elements of an operator small numberofcombinationsoftheformal seriesthatmay algebra C(x)[D], then there is a one-to-one correspondence ] give rise tohyperexponentialsolutions. C between the hyperexponential solutions y of an operator P and itsfirst orderright hand factors. In otherwords, if y is S Categories andSubject Descriptors a hyperexponential term with (q D+q ) y = 0, then y is . 1 0 · s I.1.2 [Computing Methodologies]: Symbolic and Alge- a solution of P if and only if there exist rational functions c [ braic Manipulation—Algorithms u0,...,ur−1 such that 1 General Terms P =(ur−1Dr−1+ur−2Dr−2+···+u0)(q1D+q0). v 6 Algorithms Algorithms for finding the hyperexponential solutions of a 8 linear differential equation (or equivalently, the first order 4 Keywords righthandfactorsofthecorrespondingoperators)areknown 2 sincelong. Theyareneededassubroutineinalgorithms for . Closed form solutions, D-finiteequations, Effectiveanalytic factoring operators or for finding Liouvillean solutions. See 1 continuation Chapter 4 of [11]for details and references. 0 Classicalalgorithmsfirstcompute“localsolutions”atsin- 3 1. INTRODUCTION gular points (cf. Section 2.3 below) and then test for each 1 combination of local solutions whether it gives rise to a hy- : Weconsider linear differential operators v perexponential solution. This leads to a combinatorial ex- Xi P =prDr+pr−1Dr−1+···+p0 plosionwithexponentialruntime. Thesituationissimilarto classicalalgorithmsforfactoringpolynomialsover ,which r wherep0,...,pr arepolynomialsandD representsthestan- first compute the irreducible factors modulo a prQime and a dard derivation d . Such operators act in a natural way dx then test for each combination whether it gives rise to a on elements of a differential ring containing the polynomi- factor in [x]. als. An object y is called a solution of the operator if P Q The algorithm of van Hoeij [12] avoids the combinatorial appliedtoy yieldszero. Weareinterestedinfindingthehy- explosion as follows. It picks one local solution and consid- perexponentialsolutionsofagivenoperator. Anobjecty is ers the operator Q = q D+q with q ,q C((x)) which 1 0 1 0 ∈ annihilates it. This operator is a right factor of P, though ∗Supported by the Austrian Science Fund (FWF) grant Y464-N18. notwithrationalcoefficients. Thealgorithmthenconstructs †UMR 5668 CNRS – ENS Lyon – Inria– UCBL (if possible) a left multipleB of Q with rational coefficients of order at most r 1. This leads to a nontrivial factor- − ization P = AB in C(x)[D]. The procedure is then applied recursively to A and B until a complete factorization is found. The first order factors in this factorization give Permission tomake digital orhardcopies ofall orpartofthis workfor risetoatmostr hyperexponentialcandidatesolutions(pos- personalorclassroomuseisgrantedwithoutfeeprovidedthatcopiesare siblyuptomultiplicationbyarationalfunction). Theseare notmadeordistributedforprofitorcommercialadvantageandthatcopies then checked in a second step. Van Hoeij’s algorithm re- bearthisnoticeandthefullcitationonthefirstpage.Tocopyotherwise,to minds of thepolynomial factorization algorithm of Lenstra, republish,topostonserversortoredistributetolists,requirespriorspecific Lenstra,Lovász[6,15],whichpicksonemodularfactor and permissionand/orafee. Copyright20XXACMX-XXXXX-XX-X/XX/XX...$10.00. constructs(ifpossible) amultipleofthisfactorwithinteger coefficients butsmaller degreethantheoriginal polynomial. application in the sense that we have (PQ) y =P (Q y) · · · This multiple is then a proper divisor in [x]. for all P,Q K[D] and all y E. Q ∈ ∈ The algorithm we propose below avoids the combinato- The elements y E such that P y =0 form a C-vector ∈ · rial explosion in a different way. We start from the local space V with dimV r. By making E sufficiently large it ≤ solutionsandregardthemasasymptoticexpansionsofcom- can always beassumed that dimV =r. plex functions. By means of effective analytic continuation 2.2 Hyperexponential Terms and arbitrary-precision numerical evaluation, we compute the values of these functions at some common ordinary ref- Let E be an extension of K. An element h E 0 is ∈ \{ } erence point. Then a linear algebra algorithm is used to calledhyperexponential overKifD(h)/h K. Equivalently, ∈ determineasmalllistofpossiblecombinationsoflocalsolu- hishyperexponentialifQ h=0forsomenonzerofirstorder · tions that may give rise to hyperexponential ones, possibly operator Q K[D]. ∈ up to multiplication by a rational function. These are then Two hyperexponential terms h1,h2 are called equivalent checked in a second step. Our approach was motivated by if h1/h2 K. For example, the terms exp(3x2 x) and van Hoeij’s polynomial factorization algorithm [14], which (1 2x)2e∈xp(3x2 x) areequivalent,butexp(3x2− x)and − − − associatestoeverymodularfactoracertainvectorandthen (1 2x)√2exp(3x2 x) are not. (Here and below, we use uses lattice reduction to determine a small list of combina- sta−ndard calculus n−otation to refer to elements of some ex- tions that may give rise toproper factors. tension E on which the derivation acts as the notation sug- Although our algorithm avoids the combinatorial explo- gests, e.g. D(exp(3x2 x))=(6x 1)exp(3x2 x).) sion problem, we do not claim that it runs in polynomial Every hyperexponen−tial term ca−n be written−in the form time. Indeed,nopolynomialtimealgorithmcanbeexpected h = exp( v), where v is a rational function. The additive becausethereareoperatorsP whichhavehyperexponential constantoftheintegralamountstoamultiplicativeconstant solutions y that are exponentially larger than P. Also van for h, whRich is irrelevant in our context, because P h = 0 Hoeij [12]makesnoformal statement about thecomplexity if and only if P (ch) =0 for every c C 0 . If w·e con- of his algorithm. It is clear though that his algorithm is su- sider the partial·fraction decompositio∈n of\v{a}nd integrate periortothenaivealgorithm. Similarly,webelievethatour it termwise, we obtain something of theform algorithm has chancesto outperform van Hoeij’s algorithm, n atleastinexamplesthatarenotdeliberatelydesignedtoex- g+ γ log(p ) hibitworstcaseperformance. Thereasonispartlythatdur- i i ingthecriticalcombinationphaseweonlyworkwithfloating Xi=1 pointnumbersofmoderateprecision whilevanHoeij’s algo- with g K, γ1,...,γn C and monic square free pairwise ∈ ∈ rithmingeneralneedstodoarithmeticinalgebraicnumber coprime polynomials pi C[x]. In terms of this representa- ∈ fieldswhosedegreesmaygrowduringthecomputation. An- tion,twohyperexponentialtermsareequivalentifthediffer- other advantage of our algorithm is that it is conceptually ence of the corresponding rational functions g is a constant simplerthanvanHoeij’s, atleastifwetakeforgrantedthat and any two corresponding coefficients γi differ by an inte- we can compute high-precision evaluations of D-finite func- ger. tions. The equivalence class of a hyperexponential term h is calledtheexponential part ofh. Themotivationforthister- 2. PRELIMINARIES minologyisthatwhenwearesearchingforsomehyperexpo- nentialsolution hofP andwealreadyknowitsequivalence In this section, we recall some results from the literature class, then we can take an arbitrary element h from this 0 and introduce notation that will be used in subsequent sec- class and make an ansatz h = uh for some rational func- 0 tions. tion u K. The operator P˜ :=P D D(1/h0) K[D] ∈ ⊗ − 1/h0 ∈ 2.1 DifferentialFieldsandOperatorAlgebras thenhasthepropertythatuisasolutionofP˜ ifandonlyif (cid:0) (cid:1) uh is a solution of P. This reduces the problem to finding A differential ring/field is a pair (K,D) where K is a 0 rational solutions, which is well understood and will not be ring/field and D: K K is a derivation on K, i.e., a map → discussed here [1,11]. satisfyingD(a+b)=D(a)+D(b)andD(ab)=D(a)b+aD(b) foralla,b K. Throughoutthispaper,weconsiderthedif- 2.3 Local Solutions ∈ ferential field K = C(x), where C is some (computable) Consider an operator P C(x)[D] of order r. By clear- subfield of , together with the derivation D: K K de- ∈ C → ing denominators, if necessary, we may assume that P finedbyD(c)=0forallc∈C andD(x)=1. Forsimplicity, C[x][D], say P = p Dr + +p with p = 0. A poin∈t we assume throughout that C is algebraically closed. r ··· 0 r 6 z iscalled singular ifz isarootofp ,orz= . A differential ring/field E is called an extension of K if ∈C∪{∞} r ∞ A point which is not singular is called ordinary. Note that K E,andthederivationofE restricted toK agrees with ⊆ there are only finitely many singular points, and that we the derivation of K. include the “point at infinity” always among the singular By K[D] we denote the set of all polynomials in the in- points. determinate D with coefficients in K. Addition in K[D] is Ifz=0isanordinarypointthenP admitsrlinearlyinde- defined in the usual way, and multiplication is defined sub- pendent power series solutions. If z =0 is a singular point, ject to the commutation rule Da = aD+D(a) for a K. ∈ it is still possible to find r linearly independent generalized The elements of K[D] are called operators, and they act series solutions of the form on the elements of some extension E of K in the obvious way: If P = p +p D+ +p Dr is an operator of or- m der r and y ∈0E, th1en P···y· :=r ri=1piDi(y) ∈ E. The xαexp(u(x−1/s)) bk(x1/s)log(x)k (1) noncommutativemultiplication iscompatiblewithoperator k=0 X P where α C, u C[x] with u(0) = 0, s , m analyticmeaningtoExp(e)=exp( e)=exp(u+αlogx)= ∈ ∈ ∈ N ∈ N x saonldutibo0n,.s..a,tbm0. ∈ThCe[[cxo]]m.pWuteatcioanll otfhessuechsosluoltuiotniosnsthies lwocealll- xmαaekxinpg(ua)c(hfooricesufiotrabalebrαan∈chCofatnhdeRulo∈gaCrit[xh−m1.])Eavmeroyucnhtositcoe known and will not be discussed here [13, 11]. gives rise to the same function up to some multiplicative Two series as in (1) are called equivalent if they havethe constant. same u and s and the difference of the respective values of Since Exp(e)b is a solution of P iff b is a solution of the α isin 1 . The equivalenceclasses of generalized series un- operator P (D+ e), we may assume that e = 0. Then derthisseZquivalencerelationarecalledtheexponential parts the problem⊗remainxs that the formal power series yˆ = b of the series. Adopting van Hoeij’s notation and defining may not beconvergent if 0is asingular point. However,by Exp(e):=exp( xe) for e∈C[x−1/s], we have that Exp(e1) resummation theory [2, 3] it is still possible to associate to and Exp(e ) are equivalent iff e e 1 . Note that if yân analytic function y defined on some sector 2 R 1− 2 ∈ sZ m=0and s=1,twoseries areequivalentiff theirquotient ∆=∆(d,ϕ,ρ):= z :0< z ρ d argz ϕ/2 canbeidentifiedwithaformalLaurentseries. Wewillfrom { ∈C | |≤ ∧| − |≤ } now on makenonotational distinction between Exp(e)and (with d [0,2π], ρ,ϕ > 0) such that yˆ is the asymptotic its equivalence class. expansio∈n of y for z 0 in ∆. A point z = 0 can be moved to the origin by the change Thepreciseformul→ation ofthisresultistechnicalandnot 6 of variables x˜=x z (if z C) or x˜ =1/x (if z = ). If reallyneededforourpurpose(see[3,Chap.6,10,and11]or P˜ is the operator o−btained f∈rom P by replacing x by∞x˜+z [2,Chap.5–7]forfulldetails). Itwillbemorethansufficient or 1/x˜, then a local solution of P C[x][D] at z is defined to knowthe following facts: as thelocal solution of P˜ C[x˜][D∈] at 0. Throughout the rest of∈this paper, we will use the follow- For every k = (k1,...,kq) q with k1 > >kq ing notation. P is some operator in C[x][D] of order r, by • and every d=(d1,...,dq) ∈[0,Q2π]q such that··· ∈ Wz1e,.w..r,itzen−x˜1i∈=Cxwezdien(iot=eit1s,.fi.n.i,tnesin1g)ulaanrdpox˜innts=,z1n/x=f∞or. |dj+1−dj|≤(kj−+11−kj−1)π2, j =1,...,q−1, − − thevariableswithrespecttowhichthesingularitiesatz ap- oneconstructs[3,§10.2]adifferentialsubring x i k,d C{ } pear at the origin. For i=1,...,n, we consider the vector of [[x]] [3, Theorems 51 and 53] which contains the C space Vi generated by all local solutions at zi. There may ring x of all convergent power series. C{ } besolutionswithdifferentexponentialparts,sayℓ different i parts Exp(ei,1),...,Exp(ei,ℓi) for ei,j ∈C[x˜i−1/si,j]. By • Torheemres i5s1aanddiffe5r3e]ntika,ldrifnrogmhomoxmko,rdphtiosmthe[3,geTrmhes- S C{ } V =V Exp(e )C((x˜1/si,j))[logx˜ ] of analytic functions defined on sectors of the form i,j i∩ i,j i i ∆(d ,ϕ,ρ)forsuitableϕ,ρ>0,withthepropertythat 1 we denote the vector space of all local solutions of P at zi foreveryyˆ x k,d thefunction k,d(yˆ)hasyâsits wariethwerixtpteonneVnetii,ajl(Pp)aritn(veaqnuiHvaoleeinj’tstpoa)pEerxsp[(1e3i,,ℓ1i)2.].Our Vi,j Tashyempkt,odtimc∈aepxCpc{aonn}svieorngefonrtzfo→rma0l[p3,oSw§1e0r.2se,rEiexsetrocictehe2ir]. The condition in the definition of equivalence that the sumSin theusual sense [3, Lemmas 8 and 20]. differenceofcorrespondingvaluesofαbeaninteger(rather than, say, requiring exactly the same value of α) ensures For a given operator P C[x][D] of order r, one • ∈ that the Vi,j are indeed vector spaces, because if some Vi,j can compute a tuple k and finite subsets D1,...,Dq contains, for example, thetwo series of[0,2π]suchthatanyyˆ [[x]]withP yˆ=0belongs ∈C · to x for all d as above with d / ,...,d / xα(1+x+x2+···) and xα(1+x+3x2+···) q.C{Ad}dk,idtionally, given such a d, o1ne∈cDan1 compqut∈e D then it must also contain their difference xα(2x2+ ) = ϕ,ρ>0suchthateachSk,d(yˆ)isdefinedon∆(d1,ϕ,ρ). xα+2(2+ ). ··· ··· Furthermore,givenapointz ∆(d1,ϕ,ρ),aprecision • ∈ 2.4 AnalyticSolutions ε>0,andyˆ [[x]]withP yˆ=0,onecanefficiently ∈C · compute an approximation Y of the vector Y(z) = ε atIatnisocrdlaisnsaicryalptohiantttzheformacatlupalolwyecronsevreiregsesoinluatinoenisghyˆboofuPr- (DjSk,d(yˆ))rj=−01 such that kY(z)−Yεk≤ε. ∈C hood ofz andthusgiverisetoanalyticfunctionsolutionsy The computational part of the last two items is a special of P. The correspondence is one-to-one. For any other or- caseofTheorem7ofvanderHoeven[10]. Asanapplication, dinary point z′ and a path z ; z′ avoiding singular vanderHoeven[9]showshowtofactordifferentialoperators points there exis∈tsCa matrix Mz;z′ ∈Cr×r such that usingnumerical evaluation. Notethat ourkj correspond to 1/k invanderHoeven’sarticles,andthecomponentsofthe Djy(z′) rj=−01=Mz;z′ Djy(z) rj=−01 tupjles k and d appear in reverse order. Also observethat in thelast item, z is an ordinary point, foreverysol(cid:0)utionya(cid:1)nalyticnearz.(cid:0)Therea(cid:1)realgorithms[4, sothatfromtherewecanuseeffectiveanalyticcontinuation 8] for efficiently computing the entries of Mz;z′ for any tocomputevaluesof (yˆ)anditsderivativesatanyother given polygon path z ; z′ with vertices in Q¯ to any de- ordinary point. Sk,d sired precision. In other words, we can compute arbitrary precision approximations of y and its derivatives at every 3. OUTLINE OFTHEALGORITHM ordinary point (“effectiveanalytic continuation”). Assume now that 0 is a singular point, and consider the A hyperexponential term h can be expanded as a gener- case s = 1 and m = 0, i.e., let yˆ = Exp(e)b for some e alized series at every point z , in particular at C[x−1]and b C[[x]] beaformal solution of P. To givea∈n its singularities. The resulting∈gCen∪er{a∞liz}ed series are local ∈ solutions of P if h is a solution of P. If h = exp( v) is a to the space Vi′,j′ (which possibly also has higher dimen- hyperexponential solution where v (x), and if we write sion),then,insomesense,hmustbelongtotheintersection the partial fraction decomposition o∈f vCin theformR of thevectorspaces Vi,j and Vi′,j′. e e e Ouralgorithm isbased on testing which intersections are v= 1 + 2 + + n , nontrivial. Tomaketheseintersectionsmeaningful,wemust x z1 x z2 ··· 1/x − − first map thevectorspaces wewant tointersect intoacom- wheretheei are polynomials in x˜−i 1, thenexpandingthis h mon ambient spaceW. Let E besome differentialring con- atzi yieldsageneralizedseriesinx˜i whoseexponentialpart taining C(x) as well as all the hyperexponential solutions matches Exp(ei). The components ei in the decomposition of P, and let W E be the C-vector space generated by ⊆ ofvmusthenceshowupamongtheexponentialpartsofthe solutions of P in E. Foreach i, let π besome vectorspace i local solutions of P. homomorphism IfExp(e ),...,Exp(e )are(representativesof)thedif- ferent expio,n1ential partsi,tℓihat appear among the local so- ℓi Exp(e ) ((x˜1/si,j))[logx˜ ] V πi W lutions at zi, then any hyperexponential solution must be i,j C i i ⊇ i−→ equivalent to the term exp( (e1,j1 + + en,jn)) for some Mj=1 x˜1 ··· xñ tuple(j1,...,jn). Itthenremainstocheckforeachofthese with the following properties: R candidates whether some element of its equivalence class solves the given equation. The basic structure of the algo- 1. The sum πi(Vi,1)+···+πi(Vi,ℓi) is direct. rithm for finding hyperexponential solutions is thus as follows. 2. If h ∈ W is hyperexponential, then πi−1(h) contains theformalseriesexpansionhˆ ofhatz ,possiblyupto i Algorithm 1. Input: a linear differential operator P = a multiplicative constant. p +p D+ +p Dr, p =0, with coefficients in C[x]. 0 1 r r ··· 6 Define W :=π (V ). If h is some hyperexponential solu- Output: all the hyperexponential terms h with P h=0. i,j i i,j · tion of P, say with exponentialpart 1. Let z ,...,z be the roots of p in , and let 1 n 1 r − ∈ C C e e e zn = . exp 1,j1 + 2,j2 + + n,jn , 2. For i∞=1,...,n do (cid:16)Z (cid:16) x˜1 x˜2 ··· xñ (cid:17)(cid:17) then h W for all i, and hencethe vector space W 3. Find the exponential parts Exp(ei,1),...,Exp(ei,ℓi) W∈ ii,sjinot the zero subspace (because it contai1n,js1a∩t of the local solutions of P at zi. ···∩ n,jn least h). Ourmain observationisthattherecan beatmost 4. Determine a set U 1,...,ℓ1 1,...,ℓn r tuples j =(j ,...,j ) for which W = 0 , and that they ⊆ { }×···×{ } 1 n j s.t. forevery hyperexponential solutionhequivalent to 6 { } can becomputed efficiently once we havebases of theW . exp n ei,ji we have (j ,...,j ) U. i,j i=1 x˜i 1 n ∈ Postponing the discussion of making the πi constructive 5. For(cid:0)eRacPh (j1,...,(cid:1)jn) U do tothenextsection,assumefor themoment thatW issome ∈ vector space over C, let r = dimW < be its dimension, 6. Let h0 :=exp ni=1 eix˜,iji , and compute the oper- andsupposewearegivenndifferentde∞compositionsofsub- ator P˜ :=P (D D(1/h0)). spaces of W into direct sums: ⊗(cid:0)R P− 1/h0 (cid:1) 7. C.sp.oa.mc,epuuotfehala.lrbaatsiiosn{aul1s,o.lu..ti,ounms}of⊆P˜C,(axn)doofutthpeutvuec1tho0r, WW21,,11⊕⊕WW21,,22⊕⊕······⊕⊕WW21,,ℓℓ21 ⊆⊆WW,, m 0 . . . Thereissomefreedominstep4ofthisalgorithm. Anaive approachwouldsimplybetotakeallpossiblecombinations, Wn,1⊕Wn,2⊕···⊕Wn,ℓn ⊆W. i.e., U = 1,...,ℓ 1,...,ℓ . This is a finite { 1}×··· ×{ n} Without loss of generality, we may make the following as- set, but its size is in general exponential in the number of sumptions: singular points. For finding a smaller set U, Cluzeau and van Hoeij [5] use modular techniques to quickly discard un- Each direct sum ℓi W is in fact equal to W. If necessary tuples. Ouralgorithm, explained in thefollowing • i=1 i,j not, add one more vector space to thesum. section, addresses thesameissue. ItcomputesasetU of at L most r tuples. ℓ = ℓ = = ℓ =: ℓ. If not, pad the sum with 1 2 n • ··· several copies of 0 . { } 4. THECOMBINATIONPHASE ℓ r. If not, then because the sums are supposed In general, the differential operator P may have several • ≤ tobedirect,eachdecompositionmustcontainatleast different solutions with the same exponential part, i.e., the ℓ r copies of 0 , which can bedropped. dimension of the vector spaces V might be greater than − { } i,j one. In this case, it might be that V contains some series Lemma 2. There are at most dimW =r different tuples i,j whichistheexpansionofahyperexponentialsolutionhatz i j =(j ,...,j ) 1,...,ℓ n as well as some other series which are not. If we compute 1 n ∈{ } ssioomneobfahs.isInofstVeia,jd,,weaecchanbnaositseexlepmecetnittwtoillcoinntgaeinnetrhael bexeptahne- such that Wj :=W1,j1 ∩W2,j2 ∩···∩Wn,jn 6={0}. linearcombinationofthisseriesandsomeotherone. Now,if Proof. Induction on n. For n = 1, there are only ℓ r ≤ the expansion of h at some other singular point zi′ belongs different tuples altogether: (1),(2),...,(ℓ), so the claim is obviously true. Suppose now that the claim is shown for Similarly,ifj =j ,thenW W = 0 bytheassump- 1 6 2 i,j1∩ i,j2 { } the case when n 1 decompositions of some vector space tion that W = ℓ W is a direct sum. Therefore are given. Let U− 1,...,ℓ n be a set of tuples j with j=1 i,j ⊂ { } Wj 6= {0}. Partition the elements of U according to their Wk1 ∩WLk2 =(Wu1∩Wi,j1)∩(Wu2 ∩Wi,j2) first components, . . . =Wu1 ∩Wu2∩{0}={0}. U =U1 ∪U2 ∪···∪Uℓ, This completes the proof of the loop invariant k1 = k2 6 ⇒ i.e., Uk is the set of all tuples j whose first component is k, Wk1 ∩Wk2 ={0}. for k=1,...,ℓ. Aconsequenceofthisinvariantisthat k UdimWk ≤r (Tpilhe=Feorf1ore,frao.lr.tleh.j,,en(=jm−2(,ok.1d.,,i.jfij2,e,j=dn..)p1.∈,r,o.j{.nb.1)le,,∈mℓ.).U.iwn,kℓit}pwhnlea−cWh1eaii′vo,sjefaW:{=v0ia},ljWi=6d(ii+Ws=o1,ljju1⊆t∩,i.oW.Wn.1,1t,n,ukk-,. wbiinneeteecavroslemserocyptuhiiottaneevrdeaotfuPiotswinnℓj.o=g1SsnudinobimcsmepWoatrhcieee,jstshou≤afmnWrWionf=edviLmePreyℓjn=si1∈tioWenrais,tjdio1ins,.dd2irTceahcnte, j =1,...,ℓ). By induction hypothesis, since the Wi′,j form n−1 decompositions of the space W1,k, there are at most min(r,d1+d2)2max(r,d1+d2) dimW tuples (j ,...,j ) with W = 0 . Con- sequen1t,lky, there ar2e altognether at m(oj2s,t...,jmℓ) 6 di{m}W = operations in C. For the total cost of the algorithm we dimW =r different tuplesfor theoriginal skp=a1ce W. 1,k thereforeobtain,writingUi forthesetU intheithiteration and d :=dimW and d :=dimW , P k k i,j i,j Thedesiredindextuplescanbecomputedefficientlyusing n ℓ dynamicprogramming, asshowninthefollowing algorithm. min(r,d +d )2max(r,d +d ) k i,j k i,j Algorithm3. Input: avectorspaceW ofdimensionr,and Xi=2Xj=1kX∈Ui ≤(dk+di,j)2 ≤2r a collection of subspaces W (i = 1,...,n; j = 1,...,ℓ) Osuucthputth:atthWe =set Uℓj=o1fWail,ljtfuoip,rjleis=j1,=...(,jn,.a.n.d,jℓ≤) wr.ith the ≤2r n ℓ | d2k+{2zdkdi,j+}|d2i,j {z } property W = Ln W = 0 . 1 n Xi=2Xj=1kX∈Ui(cid:0) (cid:1) j i=1 i,ji 6 { } n ℓ 1. U :={(j)T:W1,j 6={0}} ≤2r r2+2r2di,j+rd2i,j 2. For i=2,...,n do Xi=2Xj=1(cid:0) (cid:1) n 3. Unew := 2r ℓr2+2r3+r3 ∅ ≤ 4. For j =1,...,ℓ do i=2 X(cid:0) (cid:1) 5. For k U do 8nr4. ≤ ∈ 6. If W W = 0 then In the second step, we have used the bounds d2 7. Unkew∩:=iU,jn6ew{∪}{append(k,j)} Lr2emamnda2|U,ir|es≤pecrt,ivewlyh.icIhntfohlelotwhirfdrosmtep,wke∈UuisedPdktkh≤∈eUibrokuann≤dd 8. U :=Unew ℓ d2 r2, which follows from ℓPd r. 9. Return U j=1 i,j ≤ j=1 i,j ≤ PIftheobjectiveisjusttoshowthaPt thealgorithm runsin Theorem4. Algorithm3iscorrectandneedsnomorethan polynomial time, a simpler argument applies. It suffices to 8nr4 operations in C, if the bases of the Wk are cached. observethatalltheintersectionscanbedonewithanumber of operations which is at most cubic in r, then taking also Proof. Correctness is obvious by line 6 and the fact that into account that we always have U r by Lemma 2, the whenever k = (k1,...,kn) is such that Wk 6= {0} then we bound O(nℓr4)=O(nr5) follows im| m|≤ediately. necessarily also haveW = 0 . (k1,...,kn−1) 6 { } For the complexity, we first show that it is a loop invari- 5. NUMERICAL EVALUATIONATA ant that W W = 0 for any two distinct k ,k U. k1 ∩ k2 { } 1 2 ∈ This is clear for i = 1 by line 1 and the assumption in REFERENCE POINT thealgorithm specification thatW = ℓj=1W1,j is adirect Wenowturntothequestionofhowtoconstructthemor- sum. Assume it is true for some i and consider the situa- phisms π . The basic idea is to choose a reference point z L i 0 tion right before line 8. At this point, for any two distinct that is an ordinary point of P, and let W be the space of tuples k1,k2 ∈ U we have Wk1 ∩Wk2 = {0} by induction analytic solutions of theequation in a neighborhood of z0. hypothesis. We have to show that the same is true for any For each singular point z , let ∆ be a i i ptwleosdhiasvtienctthetufpolrems kk11,k=2 ∈(uU1,nje1w)., kB2y=lin(eu72,,ja2n)yfosurcshomtue- speocwtoerr sreoroietesdapapteazirinfogrinwhthicehgaenllerfaolrimzeadl z3 z0 udi1s,tuin2ct∈ifUu1a6=nduj21o,rj2j1∈6={j12,...I.f,uℓ1}.6=Tuh2e,tthuepnlebsykin1d,ku2ctaiorne sperreiteastisoonlutaisonasnaolfyPticatfuznictaidomnsitvaian isnotmere- z1 z2 hypothesis Wu1∩Wu2 ={0}, and therefore also operator Sk,d (depending on i, but not z on the series), as described in Section 2.4. 4 Wk1 ∩Wk2 =(Wu1 ∩Wi,j1)∩(Wu2 ∩Wi,j2) Such sectors exist and can be computed = 0 W W = 0 . { }∩ i,j1 ∩ i,j2 { } explicitly. Next, let γi (i = 1,...,n) be polygonal paths from z to z avoidingsingular points andleaving thestart- 2. Let h W be hyperexponential. Then the expansion i 0 pointthrough∆ (meaningthatforsomeε>0allthepoints hˆ ofhatz ∈isclearly alocal solution,sohˆ V forsomej. i i i,j on γi with adistancetozi lessthanεshould belongto∆i). We show that πi(hˆ) = ch for some c ∈. The map πi is Such paths exist. The analytic interpretations of the gener- a differential homomorphism because∈ C is (as remarked k,d alizedseriessolutionsatthesingularpointszi definedin∆i in Section 2.4) andthe(formal) exponeSntialpartsExp(ei,j) admit a unique analytic continuation along the paths γi to are mapped to analytic functions satisfying thesame differ- the neighborhood of z0. ential equations. Since h is hyperexponential, it satisfies a Wedefineπi: Vi→W asfollows. LetVi0,j bethesubspace first order linear differential equation. Since hˆ is theexpan- ofVi,jconsistingofgeneralizedseries(1)withs=1andm= sion of h, it satisfies the same equations as h. Since πi is 0, and let Vi′,j be a linear complement of Vi0,j in Vi,j. If yˆ∈ a differential homomorphism, πi(hˆ) satisfies the same equa- Vi0,j i.e., ifyˆ=Exp(ei,j)bwithei,j ∈C[x˜−i 1]andb∈C[[x˜i]], tions as ˆh. Hence πi(hˆ) and h satisfy the same first-order define πi(yˆ) to be the unique analytic continuation of the differential equation. The claim follows. functionE(e ) (b)alongγ toz ,whereE(e )refersto i,j k,d i 0 i,j S z the function z 7→ exp( z0ei,j/x˜i) with some arbitrary but The definition of the maps π1,...,πm as outlined above fixed choice of the branch of the logarithm, and is as reliesonanalyticcontinuation,aconceptwhichisonlyavail- k,d describedinSection2.4R. Setπi(yˆ)=0foryˆ∈Vi′,j,Sandthen ableif C =C. Foractualcomputations, wemust work ina extend π to V by linearity. The precise values of π (V ) computablecoefficientdomain. Atthispoint,weusenumer- i i i i,j dependonthechoiceof∆ andd(whichisarbitrary,within ical approximations. By van der Hoeven’s result quoted in i thelimitsindicatedinSection2.4),but,asshownbelow,the Section 2.4,we are able tocomputefor every given yˆ Vi,j properties of these spaces used in the algorithm donot. and every given ε>0 a vectorYε (i)r with ∈ ∈Q Proposition 5. The functions πi defined above satisfy the Dkπi(yˆ)(z0) r−1 Yε <ε. two requirements imposed in Section 4: (1) π (V )+ + k=0− i i,1 ··· ∞ tπhie(nVi,πℓii−)1i(sh)acdoinretcatinssumth;e(f2o)rmifahl siseraiehsyepxepraexnpsoionnenotfiahl taetrmzi,, Ugosriintghmthe3searae(cid:13)(cid:13)ptp(cid:0)hreonxipmeartfoiornmse,(cid:1)dthweitlihnebaa(cid:13)(cid:13)rllaalgreitbhrmaeptaicrttsookfeAepl- possibly up to a multiplicative constant. track of accumulating errors during the calculations. The test in line 6 of this algorithm requires to check whether a Proof. 1. Without loss of generality, we assume z =0. Let i certain matrix has full rank. There are two possible out- yˆj ∈ Vi,j (j = 1,...,ℓi) and consider yˆ = ℓji=1yˆj. Write comes: If during the Gaussian elimination we can find in yˆj = xαjexp(uj)bj +yˆj′ where yˆj′ ∈ Vi′,j, Pthe (αj,uj) are every iteration an entry which is definitely different from pairwise distinct, u (0)=0, and b (0)=0 unless the series zero, then the rank of the matrix is definitely maximal and j j bj is zero. Writing uj = kuj,kx−k,6choose a direction θ theintersectionofthevectorspacesisdefinitelyempty. We suchthatρeiθ ∆ forsmallρandanytwounequalu1/keiθ arethenentitledtodiscardthepossibleextensionofthepar- havedifferent r∈ealiparts. P j,k tialtupleunderconsideration. Ontheotherhand,ifduring By changing x to e−iθx, we can assume that d=0. This theGaussianeliminationweencounteracolumninwhichall ttwraonpfoorlmynsoumjiianlstou cakn(ubj,ektehikeθ)sxam−ke,osnolythiaftthteheuretahlepmasretlsveosf tthheatetnhteriienstearrseecbtailolns tisharetaclloyntnaoinnemzeproty,,tohristhcaanttehiethaecrcumreaacny j j are equal. Hence,Pwe can reorder the nonzero terms in the of theapproximation was insufficient. In this case, in order expressionofyˆbyasymptoticgrowthrate,insuchawaythat to be on the safe side, we must consider the intersection as the nonzero terms come first, u = = u and Reα = nonempty and includethe corresponding tuple. 1 t 1 =Reα , while ··· Regardlessofwhichinitialaccuracyεisused,thisvariant t ··· ofAlgorithm3producesasetoftuplesthatisguaranteedto zReα1eReu1(z) zReαpeReup(z), z 0,z>0 containallcorrectones,butmaypossiblycontainadditional ≫ → ones. With sufficientlyhigh precision, thenumberof tuples for allp t+1suchthaty =0. Usingthedefinitionof π p i ≥ 6 in the output that actually have an empty intersection will and thefact that (b )(z) tends to b (0) as z 0 in the k,d j j S → drop to zero. We don’t need to know in advance which positive reals, it follows that precisionissufficientinthissense,becauseitisnotdramatic t to havesome extra tuples in the outputas long as they are z−Reα1exp( u1(z))πi(yˆ)(z)= cjbj(0)ziImαj +o(1) not too many. As a pragmatic strategy balancing precision − j=1 and output size, one might start the algorithm with some X (as z 0, z > 0) for some nonzero constants c . Since fixed precision ε and let it abort and restart with doubled j the (α→,u ) are pairwise distinct by assumption and the precision whenever U exceeds2r, say. j j | | (Reα ,u ) are equal for j =1,...,t, theImα are pairwise Observethatthenumericalapproximationisonlyusedto j j j distinct for j =1,...,t. determinethetuplesetU,andwedonotuseittosomehow Now assume that π (yˆ) = 0. Then, for all λ > 0, the reconstruct the exact symbolic hyperexponential solutions i expression tj=1cjbj(0)(λz)iImαj tends to 0 as z →0, z > fsrioonminit.tyWpiecatlhseirteufoarteiondso.n’t expect to need very high preci- 0. Choosing λ=ep for p=1,...,t,it follows that if not all P the b (0) were zero, thet t determinant j × 6. A DETAILED EXAMPLE det (epz)iImαq p,q =ziIm(α1+···+αt)det (eiImαq)p p,q Consider the operator would(cid:0)tend to zer(cid:1)o as well, which howeve(cid:0)r is not th(cid:1)e case. P =p +p D+p D2+p D3 [x][D] 0 1 2 3 Thereforebj(0)=0forj =1,...,t,andthereforeyˆj =0for ∈Q j =1,...,t, and therefore yˆ =0 for j =1,...,ℓ . where j i p0=−105x20+3570x19−58026x18+556216x17−3456830x16+ 1,2,3). The vectors yi,j(z0),Dyi,j(z0),D2yi,j(z0) to five 14810744x15−45667732x14+104614932x13−182764261x12+ decimal digits of accuracy are as follows. 249940430x11−276371642x10+257839924x9−211785148x8+ (cid:0) (cid:1) 154714472x7−95675216x6+45214304x5−13863936x4+ 200.15 70.513 156.55 1685888x3+424960x2−182784x+20480, − − − W = 322.46 , W = 46.308 , 91.322 , p1=(x−1)x(105x19−3150x18+51456x17−489796x16+ 1,1 1184.8! 1,2 −101.17! −205.47! 2938210x15−11903624x14+34247824x13−72603516x12+ h − i h − − i 116974957x11−148046826x10+153582952x9−137261696x8+ 109046080x7−75250624x6+41559168x5−16084864x4+ 30.349 12.494 77.105 3278080x3+163840x2−231424x+32768), W2,1= 48.896 , W2,2= 5.2891 , 44.216 , − p2=−4(x−2)2(x−1)3x2(30x15−693x14+7314x13−42905x12+ h 179.66 !i h 13.066! 99.931!i 155930x11−378483x10+649718x9−828795x8+820160x7− 645092x6+398200x5−182384x4+54656x3−5696x2−2944x+1024), .74285 15.580 4.5433 8p3309=x64(−x1−0927)42(xx5−+11)055x240(x145x−1064−562x538x+91+55124x922x+84−804x44−6x2756+). W3,1= −..10469169004! , W3,2= −13015..32067!, 25..69560331! , h i h i The leading coefficient p3 has 13 distinct roots in , but 30.349 2.8557 63.199 C those coming from the degree-10-factor turn out to be ap- W = 48.896 , W = .23797 , 41.308 . 4,1 4,2 parent,sowecanignorethem. Itthusremainstostudythe −179.66 ! −.57510 ! 90.353! singular points z :=0, z :=1, z :=2, and z := . h i h i 1 2 3 4 For each singular point, we find three linearly∞indepen- We now go through Algorithm 3. Start with the partial dent generalized series solutions with two distinct exponen- tuples (1) and (2) corresponding to the vector spaces W1,1 tial parts: andW1,2,respectively. Tocomputetheintersection ofW1,1 andW weapplyGaussianeliminationtothe3 2-matrix 2,1 V1,1 =Cyˆ1,1 V1,2 =Cyˆ1,2+Cyˆ1,3, whose columns are thegenerators of W1,1 and W×2,1: V = yˆ V = yˆ + yˆ , 2,1 C 2,1 2,2 C 2,2 C 2,3 200.15 30.349 200.15 30.349 V3,1 = yˆ3,1 V3,2 = yˆ3,2+ yˆ3,3, −322.46 48.896 − 0.00 C C C − −→ V = yˆ V = yˆ + yˆ 1184.8 179.66 ! 0.00 ! 4,1 C 4,1 4,2 C 4,2 C 4,3 − where The notation 0.00 refers to some complex number z with z <5 10−3,whichmayormaynotbezero,whiletheblank yˆ1,1 =exp(x1) 1− 49x+ 3372x2+ 38834x3+··· , e|n|tries·intheleft columnsignify exactzerosthathavebeen produced by the elimination. As the remaining submatrix (cid:16) (cid:17) yˆ =√x 1 x 25x3+ , doesnotcontainanyentrywhichiscertainlynonzero,were- 1,2 − − 24 ··· gardtheintersectionasnonempty,whichinthiscasemeans (cid:16) (cid:17) yˆ1,3 =√x x2− 74x3+ 392x4+··· , W1T,1he=inWte2,r1s.ecTtihoenspaWrtial tuWple (1a)nids WextendeWdto (t1u,r1n).out 1,1 2,2 1,2 2,1 (cid:16) (cid:17) ∩ ∩ to be trivial, as they have to be if we really have W = yˆ =(x 1)3+(x 1)5 4(x 1)6+ , 1,1 2,1 − − − 3 − ··· W1,2, because the sums W1,1 W1,2 and W2,1 W2,2 are ⊕ ⊕ direct. It thus remains to consider the intersection W yˆ =exp( 1 ) 1+ 1(x 1)+ 19 (x 1)3+ , 1,2∩ 2,2 x−1 2 − 120 − ··· W2,2. Applying Gaussian elimination to the 3×4-matrix yˆ =exp( 1 )(cid:16)(x 1)2+ 23(x 1)3+ , (cid:17) whosecolumnsarethegeneratorsofW1,2 andW2,2,wefind 2,3 x−1 − 30 − ··· 70.513 25.596 12.494 77.105 (cid:16) (cid:17) − − 46.308 2.1330 5.2891 44.216 yˆ3,1 =1− 43(x−2)+ 3392(x−2)2− 637834(x−2)3+··· , −101.17 5.1548 13.066 99.931! − − yˆ = 1 exp( 1 ) 1+ 11(x 2)+ , 70.513 25.596 12.494 77.105 3,2 (x−2)2 x−2 4 − ··· − −17.50 4.440 9.777 , yˆ = 1 exp( 1 )(cid:16)(x 2)3+ 1(x 2)(cid:17)4+ , −→ − 0.00 0.00 ! 3,3 (x 2)2 x 2 − 4 − ··· − − (cid:16) (cid:17) which suggests that we have W1,2 = W2,2. We extend the yˆ4,1 =x 1+3x−1+9x−2+ 739x−3+74x−4+··· , wpaerhtiaavletuUpl=e (2(1),t1o),((22,,22)). .At the end of the first iteration, { } yˆ4,2 =√(cid:0)x 1+x−1+ 32x−2+ 163x−3+··· , (cid:1) WIn the seWcond,itseorawtieone,xwteendfintdheWp(1a,1r)ti∩alWtu3,p1le=({10,}1)antdo yˆ4,2 =√x(cid:16)x3+x+ 169x−1+ 23803x−2+···(cid:17). (w1e,(11fi,,1n2)d)⊆aWnd r3e,2corWd W(1,,1s,o2)w=eWex(t1e,1n)d=(2W,21),1t.oF(u2r,t2h,e1r)maonrde (cid:16) (cid:17) record W 3,1 ⊆= W(2,2). Finally, there is a nontrivial inter- Let us choose z0 = 3 as ordinary reference point and take (2,2,1) 3,1 the branch of the logarithm for which √x is positive and section between W(2,2) and W3,2: real on the positive real axis. The example was chosen in 12.494 77.105 15.580 4.5433 such a way that all thepower series are convergent in some 5.2891 44.216 31.307 2.6503 neighborhoodoftheexpansionpoint,sothatwedonotneed 13.066 99.931 −105.26 5.9631! to worry about sectors and resummation theory but can 12.494 77.105 15.580 4.5433 use the somewhat simpler algorithm for effective analytic 27.34 89.53 1.72 continuation in the ordinary case to compute the values of −→ − − 216. 0.00 ! the analytic functions y := π (yˆ ) (i = 1,...,4; j = i,j i i,j suggests a common subspace of dimension 1 generated by Miller’s algorithm [16]numerically solvesasimilar problem, the second listed generator of W . We therefore extend butsofarwehavenotbeenabletoturntheunderlyingcon- 3,2 the partial tuple (2,2) to (2,2,2) and record W = vergence statements into explicit error bounds that would (2,2,2) [(4.5433,2.6503,5.9631)]. Attheendoftheseconditeration, yieldanalgorithmproducingoutputwithcertifiedprecision. we haveU = (1,1,2),(2,2,1),(2,2,2) . Finally, it would of course be also interesting to see an { } For the final iteration, we see by inspection that W = analogofouralgorithmforfindinghypergeometricsolutions 4,1 W =W ,soweextend(1,1,2)to(1,1,2,1). Because of linear recurrence equations with polynomial coefficients. 2,1 (1,1,2) dimW = 1 and the sums of the vector spaces are direct, Atranslationisnotimmediatebecausethereisnonotionof 4,1 the other two partial tuples cannot also have a nontrivial local solution around a finite singularity in thiscase. intersection with W , nor can W W be nontrivial. We do however4,1have W (1,1,2W) ∩ a4n,2d W 8. REFERENCES W , so the algorithm term(i2n,2a,t1e)s⊆with4,2the outp(u2t,2,U2) ⊆= [1] S.A.Abramov and K.Yu.Kvashenko.Fast algorithms 4,2 (1,1,2,1),(2,2,1,2),(2,2,2,2) . tosearch for therational solutions oflinear differential { At this point we know that e}very hyperexponential solu- equations with polynomial coefficients. In Proceedings tion of theoperator P must haveone of thefollowing three of ISSAC’91, pages 267–270, 1991. exponential parts: [2] WernerBalser. From Divergent Power Series to Analytic Functions, volume1582 of Lecture Notes in 1 1 1 exp + from (1,1,2,1) Mathematics. Springer-Verlag, 1994. (x 2)2 x x 2 − − [3] WernerBalser. Formal power series and linear (cid:16) 1 (cid:17) √xexp from (2,2,1,2) systems of meromorphic ordinary differential x 1 equations. Springer,2000. − 1 (cid:16) 1 (cid:17) √xexp + from (2,2,2,2). [4] David V.Chudnovskyand Gregory V.Chudnovsky. x 1 x 2 Computer algebra in the service of mathematical (cid:16) − − (cid:17) Following the steps of Algorithm 1, it remains to check physicsand numbertheory.In David V. Chudnovsky whethersomerational functionmultiplesofthesetermsare and Richard D. Jenks, editors, Computers in solutions of P. The important point is that we have to do Mathematics, volume125ofLecture Notes in Pure and this only for three different candidates, while the naive al- Applied Mathematics, pages 109–232, Stanford gorithmwouldhavetogothroughall24 =16combinations. University,1986. Dekker. Indeed,it turnsoutthatP hasthefollowing threehyperex- [5] Thomas Cluzeau and Mark van Hoeij. A modular ponential solutions: algorithm to computethe exponentialsolutions of a linear differential operator. Journal of Symbolic (x 1)3 1 1 1 − exp + , √xexp , Computation, 38:1043–1076, 2004. (x 2)2 x x 2 x 1 − − − [6] A.K. Lenstra, H. W. Lenstra, and L. Lovász. (x 2)x2√x(cid:16)exp 1 +(cid:17) 1 . (cid:16) (cid:17) Factoring polynomials with rational coefficients. − x 1 x 2 Annals of Mathematics, 126:515–534, 1982. (cid:16) − − (cid:17) [7] Marc Mezzarobba. NumGfun: a package for numerical 7. CONCLUDINGREMARKS and analytic computation with d-finitefunctions. In Our algorithm as described above takes advantage of the Proceedings of ISSAC’10, 2010. fact that series expansions of hyperexponential terms can- [8] Joris van derHoeven. Fast evaluation of holonomic not involve exponential terms with ramification (s > 1) or functions. Theoretical Computer Science, logarithms (m > 0), by letting the morphisms πi map all 210(1):199–216, 1999. these irrelevant series solutions to zero. As a result, we get [9] Joris van derHoeven. Aroundthe numeric-symbolic smaller vector spaces Wi,j, which not only reduces the ex- computation of differential galois groups. Journal of pected computation time per vector space intersection but Symbolic Computation, 42:236–264, 2007. also makes it somehow more likely for intersections to be [10] Joris van derHoeven. Efficient accelero-summation of empty,thusdecreasingthechancesofgettingtuplesthatdo holonomic functions. Journal of Symbolic not correspond to hyperexponentialsolutions. Computation, 42(4):389–428, 2007. Asafurtherrefinementinthisdirection,itwouldbedesir- [11] MariusvanderPutandMichaelSinger.Galois Theory abletoexploit thefact thatifˆh=Exp(e)bistheexpansion of Linear Differential Equations. Springer, 2003. of some hyperexponential term h, then the formal power [12] Mark van Hoeij. Factorization of differential operators seriesbmustbeconvergentinsomeneighborhoodoftheex- with rational functions coefficients. Journal of pansionpoint. InsteadofthevectorspacesW usedabove, i,j Symbolic Computation, 24:537–561, 1997. it would besufficient toconsiderthesubspacesWi′,j ⊆Wi,j [13] Mark van Hoeij. Formal solutions and factorization of corresponding to generalized series solutions involving only differential operators with power series coefficients. convergent power series. Besides the advantage of having Journal of Symbolic Computation, 24(1):1–30, 1997. to work with even smaller vector spaces, an additional ad- [14] Mark van Hoeij. Factoring polynomials and the vantage would be that the numerical evaluation becomes knapsack problem. Journal of Number Theory, simpler because algorithms for the regular case [4, 8] be- 95:167–189, 2002. come applicable. Implementations of these algorithms are available [7], which to ourknowledge is not yet thecase for [15] Joachim von zurGathen and Jürgen Gerhard. Modern van der Hoeven’s general algorithm for the divergent case Computer Algebra. Cambridge UniversityPress, 1999. [10]. Unfortunately however, it is not obvious how to com- [16] Jet Wimp. Computing with Recurrence Relations. putefrom agivenbasisofWi,j abasisofthesubspaceWi′,j. Pitman Publishing Ltd., 1984.